{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GLM: Robust Linear Regression\n",
    "\n",
    "This tutorial first appeard as a post in small series on Bayesian GLMs on:\n",
    "\n",
    "  1. [The Inference Button: Bayesian GLMs made easy with PyMC3](http://twiecki.github.com/blog/2013/08/12/bayesian-glms-1/)\n",
    "  2. [This world is far from Normal(ly distributed): Robust Regression in PyMC3](http://twiecki.github.io/blog/2013/08/27/bayesian-glms-2/)\n",
    "  3. [The Best Of Both Worlds: Hierarchical Linear Regression in PyMC3](http://twiecki.github.io/blog/2014/03/17/bayesian-glms-3/)\n",
    "  \n",
    "In this blog post I will write about:\n",
    "\n",
    " - How a few outliers can largely affect the fit of linear regression models.\n",
    " - How replacing the normal likelihood with Student T distribution produces robust regression.\n",
    " - How this can easily be done with `PyMC3` and its new `glm` module by passing a `family` object.\n",
    " \n",
    "This is the second part of a series on Bayesian GLMs (click [here for part I about linear regression](http://twiecki.github.io/blog/2013/08/12/bayesian-glms-1/)). In this prior post I described how minimizing the squared distance of the regression line is the same as maximizing the likelihood of a Normal distribution with the mean coming from the regression line. This latter probabilistic expression allows us to easily formulate a Bayesian linear regression model.\n",
    "\n",
    "This worked splendidly on simulated data. The problem with simulated data though is that it's, well, simulated. In the real world things tend to get more messy and assumptions like normality are easily violated by a few outliers. \n",
    "\n",
    "Lets see what happens if we add some outliers to our simulated data from the last post."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, import our modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pymc3 as pm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import theano"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create some toy data but also add some outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 100\n",
    "true_intercept = 1\n",
    "true_slope = 2\n",
    "\n",
    "x = np.linspace(0, 1, size)\n",
    "# y = a + b*x\n",
    "true_regression_line = true_intercept + true_slope * x\n",
    "# add noise\n",
    "y = true_regression_line + np.random.normal(scale=.5, size=size)\n",
    "\n",
    "# Add outliers\n",
    "x_out = np.append(x, [.1, .15, .2])\n",
    "y_out = np.append(y, [8, 6, 9])\n",
    "\n",
    "data = dict(x=x_out, y=y_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data together with the true regression line (the three points in the upper left corner are the outliers we added)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAG5CAYAAAAEdtrhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXJ+yboATZAgKRPeyBYBBBrWiV6tVqaym26nUp1p/+an9ob6/3ar3tvb1t9bbWVkTtammr1t5abV0pqET2TUAUgghhTVD2Pfn+/pjJZGayzSSznDPzfj4ePEgyZ875zsnkvOe7nO/XnHOIiIj4WU66CyAiItJcCjMREfE9hZmIiPiewkxERHxPYSYiIr6nMBMREd9TmIlvmdmNZvZOHNtvNbPPJLNMyWBmD5rZMyk61gIzu6UZz3dmdm4M233bzJ5q6nHSJZ7zE+u5kMRQmGUYM7vezJaY2REz2xv8+g4zs3SXLVpzL5zJpAtRcjnn/tM558nfvfiTwiyDmNk3gZ8APwR6AN2BrwGTgNYpLkvLVB5P0ku/b0k3hVmGMLPOwEPAHc65551zh1zAKufcl51zJ4LbtTGzH5nZNjPbY2ZzzKxd8LGpZlZmZt8M1up2mdlNYceI5bn3mdlu4JdmdqaZvWRm5Wb2afDrvOD23wMmA4+Z2WEzeyz48yFm9rqZfWJmH5jZF8KO39XMXjSzg2a2FMhv5JzcYGYfm9k+M/vXqMcmmNm7ZrY/+DofM7PWwcfeCm62Jli2Lzb0Wuo59rfMrNTMDpnZBjO7OuyxG83sneC5/NTMPjKzz4Y93t/MFgaf+zqQ28BxajW1htcqzexXZvYzM3s5uL8lZpYftu0lZrbRzA4EfwcWta+bzez9YDlfNbNzoo7zdTPbBGyKet744HukRdjPrjGzNcGvQ02nZtYvuK+vBt9bFeG/LzNrZ2a/DpbhfTO718zKGjgnzgKtEZuCr/k/zCzfzEqC751nq3/Xwe1vNbPNwffci2bWKxHnR1LMOad/GfAPuAw4DbRsZLv/AV4EzgI6AX8F/iv42NTgPh4CWgGXA0eBM+N47n8DbYB2QFfg80D74PbPAf8bVpYFwC1h33cAtgM3AS2BMUAFMCz4+B+AZ4PbFQA7gHfqeZ3DgMPABcHyPBIs32eCj48DJgaP0w94H/i/Yc93wLlh3zf4Wuo4/nVALwIfGL8IHAF6Bh+7ETgF3Aq0AGYBOwELPv5usLxtguU/BDxTz3FujD4H4WUHfgXsAyYEX+vvgD8EH8sN7vva4O/7G8FzdEvw8auAzcDQ4HPvB0qijvN68P3Qro5jbwA+G7b9n4FvBr9+sPo1Bc+/A54Mvm9GASeAocHHvw8sBM4E8oC1QFkD594BfwHOAIYH9/UmMADoHCzXV4PbXkTgPTY2eL5/CryVwPNzbn3l1L8EXwPTXQD9S9AvEmYCu6N+VgLsB44FL4oWvKjmh21zHvBR8OupwW1bhj2+l8BFP5bnngTaNlDG0cCnYd8vIDLMvgi8HfWcJ4AHCFz0TwFDwh77T+oPs38neNEOft8hWL7P1LP9/wX+HPZ9gxei6NcSw+9nNXBV8Osbgc1hj7UPHq8H0Dd4wewQ9vg8mhdmT4U9djmwMfj1V4DFYY8ZUBZ2sf478M9hj+cQ+HBzTthxLmrg2PcBvwt+fVbwudWB/iC1wywvbD9LgeuDX28BLg177BYaD7NJYd+vAO4L+/5h4MfBr58GfhD2WMfg+6xfgs6PwixF/9TOnTn2Ablm1tI5dxrAOVcMEGySyQG6EbhwrrCa8SBGIChC+6l+ftBRAn/gsTy33Dl3PPSgWXsCtbnLCHyqBuhkZi2cc5V1vIZzgCIz2x/2s5bAb4PHb0mg5lbt47pPBRCoFYW2dc4dMbN9YWUbRKD2Uxh8XS0JXPTqFO9rMbOvAPcQuChC4ByGNxfuDivb0eA5rd7mU+fckajX2aeB19qY3WFfV/8+ofY5cmYWfn7PAX5iZg+H/cyA3tSc+/Dtoz0DvG9mHYAvEPigsqu55WzkmNX2hH19rI7ve4Tte2X1A865w8H3Se/o4zbx/EiKqM8sc7xLoDnlqga2qSDwhzzcOdcl+K+zc65jA8+J57nRSzB8ExgMFDnnziBQO4Safofo7bcDC8P238U519E5NwsoJ1BjCb+o922gvLvCtw2GUdewxx8HNgIDg2X7NlH9IXG+lpBgv8mTwJ1AV+dcF2BdI/sPL/eZwQCo1tDrPEIgjKuP3aOBbes6Vvg5MiLP73bg9qjfRzvnXEnYNvUuu+Gc20HgfXkNcAOBDyVNsYtA82K15gR7tJ0EQgmA4HnvSqAJOxHnR1JEYZYhnHP7ge8APzeza82sk5nlmNloAk1sOOeqCFxk/8fMzgYws95mdmkM+2/KczsRCMD9ZnYWgebCcHsI9GNUewkYZIGBG62C/8ab2dBg7ecF4EEza29mw4CvNnDs54HpZnZ+sLP/ISLf752Ag8BhMxtCoN+qobI19lrCdSBwkS8HsMAgmoIGtg9xzn0MLAe+Y2atzex84HMNPGUNMNzMRptZWwLNd7F6OfjcaywwGvEuamosAHOAfzGz4cHX0dnMrotj/wC/Ae4FRhD4/TXFs8FynGlmvQl8SEiU3wM3Bc9fGwJN10ucc1tJzfmRBFGYZRDn3A8ING3dS+BivIdAn9N9BPrPCH69GVhsZgeBNwjUOGIR73N/TKBDvwJYDLwS9fhPgGuDI8Eedc4dAqYB1xP4xLybmgElELiIdQz+/FfAL+s7sHNuPfB1Av1Nu4BPCfR3VPt/wAwCHfxPAn+M2sWDwK8tMNrxCzG8lvBjbyDQL/Mugd/BCGBRfdvXYQZQBHxCIDR/08CxPiQQ1G8QGFEY803kzrkKAgNVvk+gmXpgeDmdc38mcP7/EPx9rwM+W8euGvJnAjWfPzvnjsb53GoPEfjdfUTgdT5PoBWi2ZxzbwD/BvyJwPskn8D7L1XnRxKkevSUiEhSmFkpgea4NxK0v1kEBodMScT+JDOoZiYiSWNmnyfQ5Dq/GfvoaWaTgs3mgwn0X/45UWWUzKDRjCKSFGa2gMD9fjcE+1ybqjWB5vL+BG41+QPw82YXUDKKmhlFRMT31MwoIiK+56lmxtzcXNevX790F0NERDxixYoVFc65bo1t56kw69evH8uXL093MURExCPMLKbZVNTMKCIivqcwExER31OYiYiI73mqz0xEssepU6coKyvj+PHjjW8sGa9t27bk5eXRqlWrJj1fYSYiaVFWVkanTp3o168fYcsKSRZyzrFv3z7Kysro379/k/ahZkYRSYvjx4/TtWtXBZlgZnTt2rVZtXSFmYikjYJMqjX3vaAwExER31OYiYikwNSpU+OaFGLBggVMnz49Ifv98Y9/zNGjTV1Ozh8UZiLieXMWllJSWhHxs5LSCuYsLE1TifxFYSYi4gEj8zpz57xVoUArKa3gznmrGJnXucn7PHLkCFdccQWjRo2ioKCAP/4xsNj4Qw89xPjx4ykoKOC2226jemWRqVOn8o1vfIPCwkKGDh3KsmXLuOaaaxg4cCD3338/AFu3bmXIkCF8+ctfZujQoVx77bV1hshrr73Geeedx9ixY7nuuus4fPgwAK+88gpDhgxh7NixvPDCC3WW+9ixY1x//fUMHTqUq6++mmPHjoUemzVrFoWFhQwfPpwHHngAgEcffZSdO3dy4YUXcuGFF9a7ne855zzzb9y4cU5EssOGDRvi2n7R5nI35qHX3MOvbnRjHnrNLdpc3qzjP//88+6WW24Jfb9//37nnHP79u0L/WzmzJnuxRdfdM45N2XKFHfvvfc655z78Y9/7Hr27Ol27tzpjh8/7nr37u0qKircRx995AD3zjvvOOecu+mmm9wPf/jD0POXLVvmysvL3eTJk93hw4edc859//vfd9/5znfcsWPHXF5envvwww9dVVWVu+6669wVV1xRq9wPP/ywu+mmm5xzzq1Zs8a1aNHCLVu2LKLsp0+fdlOmTHFr1qxxzjl3zjnnuPLymvNV33bpVtd7AljuYsgP1cxExBeK83OZWdSXR+dvZmZRX4rzc5u1vxEjRvD6669z33338fbbb9O5c6CW949//IOioiJGjBjB/PnzWb9+feg5V155Zei5w4cPp2fPnrRp04YBAwawfft2APr06cOkSZMAmDlzJu+8807EcRcvXsyGDRuYNGkSo0eP5te//jUff/wxGzdupH///gwcOBAzY+bMmXWW+6233go9NnLkSEaOHBl67Nlnn2Xs2LGMGTOG9evXs2HDhjr3Eet2fqIwy3Dqa5BMUVJawTNLtnHXRefyzJJttd7X8Ro0aBArV65kxIgR3H///Tz00EMcP36cO+64g+eff5733nuPW2+9NeLepzZt2gCQk5MT+rr6+9OnTwO1h5hHf++c45JLLmH16tWsXr2aDRs28PTTTzfrtQB89NFH/OhHP+LNN99k7dq1XHHFFXXetxXrdn6jMMtwyehrEEm16vftYzPGcM+0wTw2Y0zE+7opdu7cSfv27Zk5cyazZ89m5cqVoYt6bm4uhw8f5vnnn497v9u2bePdd98FYN68eZx//vkRj0+cOJFFixaxefNmINB39+GHHzJkyBC2bt1KaWngg+bvf//7Ovd/wQUXMG/ePADWrVvH2rVrATh48CAdOnSgc+fO7Nmzh7///e+h53Tq1IlDhw41up2faTqrDFecnxv6w59Z1JdnlmzjsRljmt1EI5JKa8sORLxvq9/Xa8sONPm9/N577zF79mxycnJo1aoVjz/+OF26dOHWW2+loKCAHj16MH78+Lj3O3jwYH72s59x8803M2zYMGbNmhXxeLdu3fjVr37Fl770JU6cOAHAd7/7XQYNGsTcuXO54ooraN++PZMnTw4FULhZs2Zx0003MXToUIYOHcq4ceMAGDVqFGPGjGHIkCERTZ0At912G5dddhm9evXiH//4R73b+Zm54EgdLygsLHRanDM5HnntAx6dv5m7LjqXe6YNTndxRHj//fcZOnRououRUFu3bmX69OmsW7cu3UXxpbreE2a2wjlX2Nhz1cyYBRLd1yAi4jUKswyXjL4GEalbv379VCtLE4VZhmuor0FEJFNoAEiG+9qU/Fo/K87P1QAQEckoqpmJiIjvKcxERMT3FGYikpX279/Pz3/+83QXI+WKi4sTsp8bb7wxdFP5LbfckvYpsRRmIpKVGgqz6qmpEiUR+6usrExASaCkpCQh+wn31FNPMWzYsITvNx4KMxHJSt/61rcoLS1l9OjRzJ49mwULFjB58mSuvPJKhg0bxtatWykoKAht/6Mf/YgHH3wQgNLSUi677DLGjRvH5MmT2bhxY639P/jgg9xwww1MmjSJG264gcrKSmbPns348eMZOXIkTzzxBABVVVXccccdDBkyhEsuuYTLL788VOPp168f9913H2PHjuW5556r97jPPfccBQUFjBo1igsuuACA9evXM2HCBEaPHs3IkSPZtGkTAB07dgQCc0TOnj2bgoICRowYEVoCZ8GCBUydOpVrr702tJxNY5NrhC8Q2rFjR/71X/+VUaNGMXHiRPbs2QNAeXk5n//85xk/fjzjx49n0aJF8f/SGqDRjCKSfg8maa7QB+u/BeX73/8+69atY/Xq1UDgIr5y5UrWrVtH//792bp1a73Pve2225gzZw4DBw5kyZIl3HHHHcyfP7/Wdhs2bOCdd96hXbt2zJ07l86dO7Ns2TJOnDjBpEmTmDZtGitWrGDr1q1s2LCBvXv3MnToUG6++ebQPrp27crKlSsBuPjii+s87kMPPcSrr75K79692b9/PwBz5szh7rvv5stf/jInT56sVbN74YUXWL16NWvWrKGiooLx48eHgnDVqlWsX7+eXr16MWnSJBYtWlRrjsn6HDlyhIkTJ/K9732Pe++9lyeffJL777+fu+++m2984xucf/75bNu2jUsvvZT3338/pn3GQmEmIhI0YcIE+vfv3+A2hw8fpqSkhOuuuy70s+o5FqNdeeWVtGvXDggsyLl27dpQrevAgQNs2rSJd955h+uuu46cnBx69OgRWkCz2he/+MVGjztp0iRuvPFGvvCFL3DNNdcAcN555/G9732PsrKy0CKi4d555x2+9KUv0aJFC7p3786UKVNYtmwZZ5xxBhMmTCAvLw+A0aNHs3Xr1pjDrHXr1kyfPh2AcePG8frrrwPwxhtvRPSrHTx4kMOHD4dqis2lMBOR9GugBpVKHTp0CH3dsmVLqqqqQt9Xz6hfVVVFly5dQjW6WPfnnOOnP/0pl156acQ2f/vb32LaR0PHnTNnDkuWLOHll19m3LhxrFixghkzZlBUVMTLL7/M5ZdfzhNPPMFFF13UaJmBiOVtWrRoEVefX6tWrULL3oQ/t6qqisWLF9O2bduY9xUP9ZmJSFYKXxalLt27d2fv3r3s27ePEydO8NJLLwFwxhln0L9/f5577jkgEFJr1qxp9HiXXnopjz/+OKdOnQLgww8/5MiRI0yaNIk//elPVFVVsWfPHhYsWFDn8xs6bmlpKUVFRTz00EN069aN7du3s2XLFgYMGMBdd93FVVddFVoqptrkyZP54x//SGVlJeXl5bz11ltMmDCh0dfRVNOmTeOnP/1p6PtYPgzEQ2EmIlmpa9euTJo0iYKCAmbPnl3r8VatWvHv//7vTJgwgUsuuYQhQ4aEHvvd737H008/zahRoxg+fDh/+ctfGj3eLbfcwrBhwxg7diwFBQXcfvvtnD59ms9//vPk5eUxbNgwZs6cydixY0OrXker77izZ89mxIgRFBQUUFxczKhRo3j22WcpKChg9OjRrFu3jq985SsR+7r66qsZOXIko0aN4qKLLuIHP/gBPXr0iOcUxuXRRx9l+fLljBw5kmHDhjFnzpyE7l9LwIhIWmTiEjBNVd13tG/fPiZMmMCiRYuSGixe1ZwlYNRnJiKSZtOnT2f//v2cPHmSf/u3f8vKIGsuhZmISJrV108msVOfmYikjZe6OSS9mvteUJiJSFq0bduWffv2KdAE5xz79u1r1rB9NTOKSFrk5eVRVlZGeXl5uosiHtC2bdvQjdpNoTATkbRo1apVo7NtiMRKzYwiIuJ7CjMREfE9hZmIiPiewkxERHxPYSYiIr6nMBMREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER30tqmJnZN8xsvZmtM7Pfm1nTV17zkTkLSykprYj4WUlpBXMWlqapRCIimS1pYWZmvYG7gELnXAHQArg+WcfzkpF5nblz3qpQoJWUVnDnvFWMzOuc5pKJiGSmZC/O2RJoZ2angPbAziQfzxOK83N5bMYY7py3iplFfXlmyTYemzGG4vzcdBdNRCQjJa1m5pzbAfwI2AbsAg44516L3s7MbjOz5Wa2PJOWTy/Oz2VmUV8enb+ZmUV9FWQiIkmUzGbGM4GrgP5AL6CDmc2M3s45N9c5V+icK+zWrVuyipNyJaUVPLNkG3dddC7PLNlWqw9NREQSJ5kDQD4DfOScK3fOnQJeAIqTeDzPqO4je2zGGO6ZNjjU5KhAExFJjmSG2TZgopm1NzMDLgbeT+LxPGNt2YGIPrLqPrS1ZQfSXDIRkcyUtAEgzrklZvY8sBI4DawC5ibreF7ytSn5tX5WnJ+rfjMRkSRJ6mhG59wDwAPJPIaIiIhmABEREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER31OYiYiI7ynMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER8T2EmIiK+pzATERHfU5iJiIjvKcxERMT3FGYiIuJ7CjMREfE9hZmIiPiewkxERHxPYSYiIr6nMBMREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER31OYiYiI7ynMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER8T2EmIiK+pzATERHfU5iJiIjvKcxERMT3FGYiIuJ7CjMREfE9hZmIiPiewkxERHxPYSYiIr6nMBMREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER31OYiYiI7ynMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER8T2EmIiK+pzATERHfU5iJiIjvKcxERMT3FGYiIuJ7CjMREfE9hZmIiPiewkxERHwvqWFmZl3M7Hkz22hm75vZeck8XqaZs7CUktKKiJ+VlFYwZ2FpmkokIuJNya6Z/QR4xTk3BBgFvJ/k42WUkXmduXPeqlCglZRWcOe8VYzM65zmkomIeEvLZO3YzDoDFwA3AjjnTgInk3W8TFScn8tjM8Zw57xVzCzqyzNLtvHYjDEU5+emu2giIp6SzJpZf6Ac+KWZrTKzp8ysQ/RGZnabmS03s+Xl5eVJLI4/FefnMrOoL4/O38zMor4KMhGROiQzzFoCY4HHnXNjgCPAt6I3cs7Ndc4VOucKu3XrlsTi+FNJaQXPLNnGXRedyzNLttXqQxMRkeSGWRlQ5pxbEvz+eQLhJjGq7iN7bMYY7pk2ONTkqEATEYmUtDBzzu0GtpvZ4OCPLgY2JOt4mWht2YGIPrLqPrS1ZQfSXDIREW8x51zydm42GngKaA1sAW5yzn1a3/aFhYVu+fLlSSuPiIj4i5mtcM4VNrZd0kYzAjjnVgONFkJERKQ5NAOIiIj4nsJMRER8T2EmIiK+pzATERHfU5iJiIjvKcxERMT3FGYiIuJ7CjMREfE9hZmIiPiewkxERHxPYSYiIr6nMBMREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER31OYiYiI7ynMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER8T2EmIiK+pzATERHfU5iJiIjvKcxERMT3FGYiIuJ7CjMREfE9hZmIiPiewkxERHxPYSYiIr6nMBMREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER31OYiYiI7ynMRETE9xoNMzP7P2Z2ZioKIyIi0hSx1My6A8vM7Fkzu8zMLNmFEhERiUejYeacux8YCDwN3AhsMrP/NLP8JJdNREQkJjH1mTnnHLA7+O80cCbwvJn9IIllExERiUnLxjYws7uBrwAVwFPAbOfcKTPLATYB9ya3iCIiIg1rNMyAs4BrnHMfh//QOVdlZtOTUywREZHYNRpmzrkHGnjs/cQWR0REJH66z0xERHxPYSYiIr6nMBMREd9TmImIiO8pzERExPcUZiIi4nsKMxER8T2FmYiI+J7CTEREfE9hJiIivqcwExER31OYiYiI7ynMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER8T2EmIiK+l/QwM7MWZrbKzF5K9rFERCQ7paJmdjfwfgqOIyIiWSqpYWZmecAVwFPJPI6IiGS3ZNfMfgzcC1TVt4GZ3WZmy81seXl5eZKLIyIimShpYWZm04G9zrkVDW3nnJvrnCt0zhV269YtWcUREZEMlsya2STgSjPbCvwBuMjMnkni8UREJEslLcycc//inMtzzvUDrgfmO+dmJut4IiKSvXSfmYiI+F7LVBzEObcAWJCKY4mISPZRzUxERHxPYSYiIr6nMBMREd9TmImIZKg5C0spKa2I+FlJaQVzFpamqUTJozATEclQI/M6c+e8VaFAKymt4M55qxiZ1znNJUu8lIxmFBGR1CvOz+WxGWO4c94qZhb15Zkl23hsxhiK83PTXbSEU81MRCSDFefnMrOoL4/O38zMor4ZGWSgMBMRyWglpRU8s2Qbd110Ls8s2VarDy1TKMxERDJUdR/ZYzPGcM+0waEmx0wMNIWZiEiGWlt2IKKPrLoPbW3ZgaQcL52jJxVmIiIZ6mtT8mv1kRXn5/K1KflJOV46R08qzERE4pBN927FK3z05COvfRBq4kzFoBOFmYhIHLLp3q2mSNfoSd1nJiISh2y6d6spokdPTszvqpqZiIgXZcu9W/FK5+hJhZmISJyy5d6teKV69GQ4c84l/SCxKiwsdMuXL093MURE6hVe+yjOz631vSSWma1wzhU2tp1qZiIicUhn7UPqp5qZiIh4lmpmIiKSNRRmIiIeoxuz46cwExHxGN2YHT/dNC0i4jG6MTt+qpmJiHiQbsyOj8JMRCQNGusXS/SN2ZneD6cwExFJg4b6xZIxLVSi+uG8GooKMxGRNGhouZRk3JidqOVZvDo4RTdNi4ik0SOvfcCj8zdz10Xncs+0wb44XnWApWJwim6aFhHxuFRPWFzf8eJtOvTi4BSFmYhkBa/19aR6uZSGjhdv06EXVw1QmIlIVvBaX0+qJyxu6Hjx9Kelc82yhqjPTEQizFlYysi8zhEXspLSCtaWHeBrU/LTWLLmly2VfT1+FEt/WqrfH+ozE5Em8VoNJlxzy+bFvh6viLXp8GtT8mudt+L83Mggq6qEPevh2KfJLHIETWclIhGaOpVSKj6xN3eap+gL9sT8rgo0ai84OjG/a8xD9+csLGXM2TkUtdoC25dC2VJOb1tKy1OH4eq5MOqLKXkNCjMRqSW8BnPXRefGdMGvrjXVtQJzussGzbtgZ7rG+tMiOAf7NsP2JbB9KTdseZd2+zcBNV1WLYETHXrRpvJEyl6D+sxEpJam9i2lok+qqcfwcl+gp504DDtXBsNrGZQtrdV8WJXTmnVV/TjVcxx/2NOLL159DYUjChJy+Fj7zFQzE5EIzanBNLXWlIqy1RVYxfm5WV8ri+AcfLoVypaFal7sWQeuKnK7jt2hTxH0mQB9isjpOYo35m8N/d4LRyT/5u9oCjMRiRBXk1OUZPdJNadsUodTx2HX6prg2r4UjuyN3MZasLfTUCp7FdJz+JRAgHXpS8mWfYFabZ98T/RFqplRRBIiutYU/b14wIEdgeAqWxYIrl1roOpU5DbtzgrWusYH/u81hpKFVmYOAAAdVElEQVTtxxrsD03m7z3WZkaFmYgkhPqkPOb0Sdj9XjC8grWugzuiNjLoPhzyxtc0G541AMxq7a6+vspk/94VZiIiCeCbkD68NzQ0nu1LYecqOH08cps2nSGvsKbm1bsQ2p4R8yFSPSkyaACIiGSBVARNqm45iEvladi7vqafq2xpYOBGtNxBkDchOFBjAuQOhpymzZXhhX6xhijMRMS3UhE0zb1ROyGOfhLq59rx3gJ6HN5Ai9NHIzY51aIdrfqOD4ZXUaAG1v6shBzeD/foKcxExLeig+bJtz/inmkDE15TS/YtBxGqqqDig8haV8WHoYd7B/8/3rEvbftPZEu74Xx7eXvuuv5Kigf2SEqR/DCKVGEmkuV80ydUj/CguXpMLx5fsIXhvTontKaW1Ca24wdhx/Kw8FoOJ6Jmzm/ZFnqNCQ3UWFZ5Lrf/eTszOwVril9Obg3JD/foKcxEspwn+4TiEB00s6YOSGiTYEKb2JyDfaXBQRrBGTX2biB8KigAzsirGRqfNwF6jICWrUMPjwdm7rTU1BR9QmEmkuWS1SeUihpffUEzZVC3hF3om9XEdvJocCqosCbDo/sit8lpBT1HRg7U6JzX4G69PhgjHRRmIpKUPqFU1PjqCppZUwfwyGubEnahj7mJzTnYvy00FdTe99+m2+EPMVcZsdnR1l1pP+C80FRQ9BwFrdrFXJ6Gaopryw74usm4OXSfmYgkbYLgVC+GmdJZSE4dD8ygEbopeRkc3h2xyWlyOH7WUDrmn8eHrYdxz7tt+PaMSyk+t1uTD9tQjbehDxB+rbnppmkRiUmyAyCVN9qGX+irvwZCNZNm1VIO7gqbCmpJIMgqT0Zu0+7MiObCxSf6ccdzHzQa5olskk3EBwgvDQrSStMiEpOG+oSaK9bVixMlfBXk6lpK9c/jWpW68hTsWAmL58DzN8P/FMAjQ+C5r8K7jwUCrfIUn3TIZ/e5X4SrfgZfXwb3fkTJxJ8zx10N/S9g4pC+Ma1sncjVvROxmraXVxuvj2pmIpIUXph4OOZaypGK4CCNYM1rx0o4fSxymzZnQO9xNXMY5hVSsuNUg68xnlpSoppkvbaf5tJ0ViIx8lKTihc19fx44UbbOge2VFUGhsNvX1rTZPjJltpP7npucGj8+EB4dRsCOS2i9k+9I0HjHdKfiEE4ibyNIKU3iieAwkyynt/vs0q2pp4fL9xoW1JawV8Wb+CRMfupePcF9m/ZSZdP3oOThyI3bNU+UOuqDq68CdCha0zHqO+iH2+YJ2K4fSI/QPht+L+aGUXwTpOKV/nm/FRVwb5NsH0Jeza8zZHNJQygrPZ2Xc6pCa0+46H7CGjRtM/2iTg3XmiS9Wp51MwoWS3epjG/NamkWqLPT8Kadk8cgh0rAsPiq/u7ju8HoHv1Ni3aQK/R0GcCG1sNZcmpc/nqtKJmlT+8zIlo1vNCk6yXyxML1cwkI8X7ydI3NY80SfT5adInf+cCfVvV/VzblwWWQXFVkdt16llzQ3LehMDsGi3bNLmsDVF/a/LpPjPJerFegL3UpJJIibrQJuv8NPr7OXkUdq5i8VuvMOjkBs76ZDUcjRran9MSeowMNhkG5zLsnFfnSsniT2pmlKwXa9OYH5tUYpGogS3JOj8Rv58L8ynuegzW/almHsPda6HqNBPDn9Q+l31njeaZHd25+JLpFIy/MK6poCRzqWYmGUtNhx4+B6dPsHb5W7z+6l/53Fnb6bxvNd35JHIby4Gzh0OfCWxqM5T/924bphRN4Jml273zOiTpVDOTrOaHlXFTwTMDWw7tjrgpuWrHKkZWnWQkQHAS+QN0pKrXOM4cPDkwwrD3OGjTCYCBwBT3Qfpfh3iWwkwyUnTT2NqyA8yaOiCiaSwbOurTcq9Q5WnYs65myZPtSwKzyYfJAY52Hkj7ARNDM2qsP3Ama3ccqvP34bd7niT1FGaSkaIviOH9R0CT+4+8pLEBHimrnR7ZVzPCsGxZYKj8qaOR27TuCHmFNZPw5hXSvt2ZEZsUd4Pic8+u9bpKSiu4/bcrmD6yJ/dMG5y1tWxpmMJMskL1oAVP9h81UfQAj395YS0vrd3FEzeMA2rXRhMxcOOJBR9S1LGc0XxQc2/XJ6W1NzxrQNhUUEVw9tBaU0HF+rr+umYnAJ8b1QvInAE6klgKM8kanuk/SpDogH5p7a6Ix6Nro9XPiet1H9sPO5aHRhj+8/ZltDx1OHITWnPy7NF0HjSpJsA6NG+UY/jrenX9Hp64YVxEuVM9LZZ4n8JMskYm9rtEB3R1E1yTap/Owb7NwRuSg8PjyzcCNSOeWwInOvTiH0f706bfRJ76uBt3fumfOG9gz6S+Lr//niT5FGaSFRLVf+S1GR/qCuiYQ+DE4UD/VvUqyWVL4dinkdu0aA09R0U0GbY5oycbwhbcTHSQ1fe6/B5oXnvvZBqFmWSFRN3466UZ9usK6Nt/uwKgdgg4B59ujRxhuKeOqaA6dq+ZCqpPUSDIoqaCSkbQhF/oq1/XrKkDqKyqWWIlkX2c6QgWL713MpFumhaJk1duRG5o1N9/fW4g7y1fyJuvvcQNebvp+skaOLI3cgfWAnoU1ARX3njo0rfBqaCSPbVV9QeMFjnw+IItEcdJZNCkawozr7x3/ERzM4ok0SNhzWz3TBuc7uLAgR28/tqLDK/8gF6H1sKutVB1KnKb9l1rljzpUwS9xkDrDnEdJpk1mlRf6NMVLPG+d7K9eVIzgIgkSTr7c+YsLGVUz3ac135XqLnwxEeLaXN0F5dEbGnQvaBmsck+RYHh8s2cgDeZC26metBH+PEm5ddeiDMZgdGU946aJ2OTtDAzsz7AbwgsK+SAuc65nyTreCKp0NSBJM36dH14b2gqqBmbS2j9j7XAydDDbYDTrTrRsm91X9d46F0Ibc9o5qtNrVR/SAg/3i9LtnL7b1eEbgFIRmA09b2TifdIJkPSmhnNrCfQ0zm30sw6ASuAf3LObajvOWpmFK9raijF3EdTeTqwRlf10PiypYGBG1E+ojdHzh7LC+W9+dz0qxgzZiLk5CTypcbE68vMxHO86sEzNxX3S0pgNPdcea5pO0U812dmZn8BHnPOvV7fNgozb8v2tvvmqrOPpmcOf3vlr4yo2kifw+/BjpVw6kjkE1t1gLxxwf6uIsgr5JF3yj1xYUtUCKX6vVXf8R6bv5mS0n1pP6/RsnngiKfCzMz6AW8BBc65g1GP3QbcBtC3b99xH3/8cdLLI02TKYtYpi2Uq6r49YuvsGHZfG7ovZuCqg+g4sNamx3v2Jc3j/RjyPjPkD/mQjh7GLSo6RHw2oXNa+VpKq++jkz5u2sqz4SZmXUEFgLfc8690NC2qpl5n1f/4OORsovD8YMRU0Gd3raUlqcORW7Tsi30GsMaG8STW7pRUHQxc1cd8d2q2H5vAvPqeQW1iHgizMysFfAS8Kpz7pHGtleYpV5T/lBSeeFK1h9ywkPZOfhkS+RUUHs3ED4VFMCJ9j051G0MT318NqsZxN0zrsa1aM2d81YxZVA3/rxqR4Pn1YsXtkz4gOPF8yoBaR+ab2YGPA28H0uQSXrEO+w31SPOkjUsudnDwE8egZ2rguEVnArq6L7IbXJaBaeCmsBrh87hrMGTKBw5gjbABaUV/O63K/jpwm1s3H2IWVMH8PiCLY2e12QOjW+KTFkE1WvnNRUyLcCTOZrxfOBt4D2ges6cbzvn/lbfc1QzS49YP1k31BSztuyAr26mjWufzgUWl6xes2v7Utj9HrjKyO06nB28p2tCYLBGrzHQqm29Zaiu4V49phcLP6zwZBNXYzLtguhHSR9hm2aeaGaMl8IsfWJpOmzoj6ahGlQi/jAS2bTZWCiP6tmW89qWBRebXMrJrYtpfaw8cifWAroPD62STN545qytZGSfLjFdVMLD9Mm3P+KeaQO5dXJ+o88TidacUPJDE7HCTGKWqDd0sv4wEr3fWqF8cBcfLH+D4x8tZsCx9bSpeI/WnI54zqnWXWh1TlFNzavXWGjTsc5yNnZR8csnYvGP5vyNeH3wTtr7zMQfEtnnkYzpiBLeJ1N5iq+dux+2vw6rggM1Dmwn/E/YYWyiL8e6j+VPe3tz1ef+ibFjxjc6FVSsMzUkagZ/SR6/NZ829W8vk5baSf2UARlqzsJSSkorIn5WUlrBnIV1LCnvIQ1dWOMV/YcRfT7SUr7D5bDxZXj9AfjFZ+G/+sCTF8Er98G6P8GB7dDmDMi/CKZ8C2b+CfvWx/y1+Hmu/PgLdD7vq4wdO6HOIKvrdw4wpEcnHp2/mZlFfesdxBH98+L83NBF0q/vpUxS3Wxe/Xuo/lA1Mq9zg89L1++uKX974R8U75k2OPRBLBF/t2nhnPPMv3Hjxjm/WrS53I156DW3aHN5nd/H6vEFm2s9Z9Hmcvf4gs0JK2syJOr1N0vlaed2rXVu6VPO/ek2534y2rkHzqj979Gxzv15lnPLfuHc7nXOVVbW+VoefnVjg6+hrtdc8MArruCBVxp9bkOScS6j31ePL9js5r61OeJ95Yf3WSrF+j6o6zmp/Dto6jH9cq0BlrsY8iPtARb+z89h5lzNm2jG3HddwQOvRLxRYn2TeCIUmiAtfxhHP3Huw9ece/O7zv3qc859r3ft4PpuD+d+cblzrz/o3Ma/O3e4osFdxnv+wy941UGWiN9dUy6kseyvej9z39rs+t33kpv71uZmlzWTPfzqRnfOfS+5h1/dGPNzEv27a4xfQqmpYg0zDQBJsOrO1LatcvjFjeOb1MHvhxFGKVdVFZj6qXqV5O3LoOKD2tt1OSdihCHdCyKmgmpMc24in5Tfla9H9Vc0p58l0R3z0e+r6nvb9D6rWyYPqvATDQBJg7qWlGjKDNypXtfJk04cgh0rImePPx7VT9aiTeBervB7uzp1b9Zh4715NrqvIp7nNiQZHfPR76tbJ+dz6Njp7H6f1aM5A48yaVCFnyjMEqSuN//Nv1rWpAtF1v0xVE8FFbopeVlgGRRXFbldp141qyTnTYCeI6Flm/SUmeTNfpHM/Ya/rzq1a5ld77M4NHXEaabMiOJHamZMkOjmqer1kUb07szG3YfibmLM6HuQTh6tmQqqbFmg5nU0agRVTkvoMbJmscm8CdClT3rKW4+GmiSBJg/tTsaw8Oj30ZNvl/KfL2/k21cM4dbJ+Zn5PksDvw3p9wPdNJ1GzQmkjPtjcA4OlNVMA1UWnAqqKvKmZNrn1vRz9Z0YnAqqXYO79vK58tqHkuhzNWdhKS1yoLKqpmnVK+dOJJzCLI1ScZH17IX89AnYtTY0FRTbl8GhnZHbWA6cPbymr6vPBDizP5jF9bq8FhjRNJBHpPk0ACSNUjEDd7Jmk4/bod3BQRrBJsOdq6HyROQ2bbvUDNDoMx56j4M2nercXTyvK9YZN9JFA3lEUkdh5mGN1VJSfiGvPAV71gVqW9U1r/21R/CROzhY4yoK/Ot6LuTENtlMvAHl5cBIxEAez9bARTxGYeZhjdVSkn4hP7Iv2FQY/LdzJZw6GrlN606QNy5Y6yoKfN3uzGYdNp7X5dWRn4ka1eaZGriIx2VEmGXqp9fGain1XcibdD6qKqF8Y01wbV8Cn9Qxn9xZ+TUDNfoUwdlDIadFQl93rAHl5WHQiZpM2OtNqSJekRFhlsmfXuurpTR0IY/pfBzbDzuW1zQZ7lgBJw5GHrxlu0D/VvXQ+D4ToEPzL6LxrIvWUEB5efb5RPaberkpVcQrMmY0Y7JHjqVraHN9r6ux2lfE8xZ/zNPTOzPGNtXUvMo3AlG/+859wvq6JgSngmqVsNcS/ZoSuWJ1ptbOQaMiJbtl3WjGZH96ja4xtMghdNMpkJTaYEO1rwY/+Z84TLGt59Heb3Dy7cXMarOFdi9G1bpatIaeoyOngjqjZ8LK3pCGms7q+r3FUqPJ1Nq5l5tSRbwkY8Is2QMB6roAf/uKITy+YAuHjp1OyifmmJrRnINPt/Lm639lWOUH9Dy4NjDi0FVxPkAL4DScbHc2rftNrKl59RyVsqmg6qo1Qc26X4n48JGpfUtebkoV8ZKMCLPmfHqNp3kquvaX7Ila66x99e1AcctNsOilmoEaR8q5OGwbZy3YwADWMIjC8y/lULdx3PqXPTxWODYtF8C6ak23/3YFQFIn0s2Ei30q7lkUyQQZEWaNfXqNZ8BBQ81TaZmo9cCOsDkMlwRm16g6FblNu7OgTxEfdyjgP9Z0orL7KJbvOsETN4xjULA8j3WsSNun+eha0y9LtgLwxA3jEtp05tVh+iKSfBkzAKQhjU17FEsHe0omaj19MjBvYWgqqKVwcEfkNpYDZw+rGRrfZwKcNQDMgOavo5TMgRTJWver+vlentpKRJom6waANKSx/pRYmqeia3+VVfDtK4ZQWRV5jLhqP4f3Rk0FtQpOH4/cpm3nQHBVD43vPQ7anlHn7hJRM0nWQIpkrftVTX1LItktK2pm1eqrtaRk6HPl6cAaXeGLTX66tfZ2uYPD7usqgtxBMU0FlciaSaLPh2pNItJUqplFqa/WkrShz0c/CVtscinsWAmnjkRu07pj8Kbk4ND4vEJof1aTDpfImkmiB1IkqmyZfC+ZiDRPVtTMom/KbZEDjy/YEvF9s25+rqqCig9qVknevgT2baq93Zn9Q/d1rWYwS4505/apgxL4ShOj+nwN7dGJtTsOhAZqVD+WrvBIRw1PASqSXrHWzGKbytznwmsGI/M68/iCLcyaOiA0mvHxBVsYmdc5tH1xfm7DF6rjB6F0Piz4Pvz2GvjvfvDzifDXu2H1M4Ega9kW+hbDpLt5f+ocPpPzFCXT34Br5lJy1j9x8yvHGdHnLOYsLKWkNHKV5ZLSCuYsrD0vYjzbNlV4QHz9onMBuP23KygprQg9Fn6uUim87/OR1z5ISVNldR9i9XlP9znwq1S8dyW7ZUXNLFpcfULOwSdbapoLty+FvRuoNRXUGXmRi012HwEtWzd6zHhqG6momUTXRKrvCRvRuzMbdx/yRD9Xc0dsxkvTSTWf+k2lqdRn1oAG+4ROHgn0b1Wvkly2FI7ui9xBTqvADBrhU0F17t2kY8Yzc0WiZrmI90bxm4r7JfVG5HhXl071vWSZeDN2qmXqDC3iHVkZZqEL4oX5vLl4OdNzFjHo5PuBWtfu98BVRj6hw9lhta6iwJyGrdo27Zh1XITjuVgm4sLanBvFkxEesZYnXfMU6mbsxNCHAkmm7AqzU8d5b/lCFr/2Eq/32UPX99ZwT9VueCdsG2sRqHVVD43vMx66nBO6KbkpGrsIx3OxTMSFNdZPyakKj1jLk457yTTRb+LoQ4EkU2b3mR3cWTPCsGwp7FoDlScjt2l3Jp+cNYZNrYdSdMFl0GsstOmYuDIQ3/pdqewza6zvKdUj+VLdFxYLjWZMDPWZSVPF2meWWWG2aw18/G7NYI2DZVEbWGBl5Op+rj4ToOu5zap1NVc8F8umXljret6Tb5fyyGubuHVyf0/0X2iQRWbThwJpquwMs99cBVsW1Hzf5ozgTclFrG85hCUn+3PzxaObXU6/Scm8kgksX7rLIyLekZ33mQ2ZDqO/TGnRd7ku5xFKrl0OX/lfSvreyg0LOjKkX166S5gW0fdnPfLaplCQhT++tuxAWsrXUF+YiEgsMqtmFiYRzVaZ1jTixT4pEZGGZGfNLEz4MOCZRX2bNYQ9E2Z/iB5JFj0bg4iIn2Xs0PxUDmH3Og0vF5FMl5E1s/CL9z3TBocCqSm1kUTU8MKlY4469UmJSKbLyDBL5MU70c1z6Wi6/NqU/Foh3OhkyiIiPpKxA0ASIVlDxnVPlYhIbLJ+AEgiJKt5LtFNlyIi2U5h1oBkNc9lw8hCrV8lIqmkMEuxRA5O8bJMuq1BRLxPfWYplmk3YjdEfYMi0lxanNOj6gqs4vzcjLzIa/0qEUkVNTNK0mRD36CIeIPCTJIiW/oGRcQbFGbSZA2NWNSsIyKSSgozabKGRixq1hERSSUNAJEmy5SJmEXE/1Qzk2bRbCYi4gUKM2kWjVgUES9QmEmTacSiiHiFwkyaTCMWRcQrNJ2ViIh4lpaAERGRrKEwExER31OYhdEaXCIi/qQwC6M1uERE/EkzgITRjBYiIv6kmlkUzWghIuI/CrMomtFCRMR/FGZhNKOFiIg/KczCaEYLERF/0gwgIiLiWZoBREREsobCTEREfE9hJiIivqcwExER30tqmJnZZWb2gZltNrNvJfNYIiKSvZIWZmbWAvgZ8FlgGPAlMxuWrOOJiEj2SmbNbAKw2Tm3xTl3EvgDcFUSjyciIlkqmWHWG9ge9n1Z8GciIiIJlfYBIGZ2m5ktN7Pl5eXl6S6OiIj4UDLDbAfQJ+z7vODPIjjn5jrnCp1zhd26dUticUREJFMlM8yWAQPNrL+ZtQauB15M4vFERCRLJW1xTufcaTO7E3gVaAH8wjm3PlnHExGR7OWpiYbNrBz4uJm7yQW0ZkvddG7qp3NTP52b+unc1C9R5+Yc51yjfVCeCrNEMLPlscywnI10buqnc1M/nZv66dzUL9XnJu2jGUVERJpLYSYiIr6XiWE2N90F8DCdm/rp3NRP56Z+Ojf1S+m5ybg+MxERyT6ZWDMTEZEsozATERHf822YNbZWmpm1MbM/Bh9fYmb9Ul/K9Ijh3NxjZhvMbK2ZvWlm56SjnOkQ6xp7ZvZ5M3NmljXDrmM5N2b2heB7Z72ZzUt1GdMlhr+pvmb2DzNbFfy7ujwd5UwHM/uFme01s3X1PG5m9mjw3K01s7FJKYhzznf/CMwoUgoMAFoDa4BhUdvcAcwJfn098Md0l9tD5+ZCoH3w61k6N7W26wS8BSwGCtNdbq+cG2AgsAo4M/j92ekut4fOzVxgVvDrYcDWdJc7hefnAmAssK6exy8H/g4YMBFYkoxy+LVmFstaaVcBvw5+/TxwsZlZCsuYLo2eG+fcP5xzR4PfLiYwCXQ2iHWNvf8A/hs4nsrCpVks5+ZW4GfOuU8BnHN7U1zGdInl3DjgjODXnYGdKSxfWjnn3gI+aWCTq4DfuIDFQBcz65nocvg1zGJZKy20jXPuNHAA6JqS0qVXvOvI/TOBT03ZoNFzE2wC6eOcezmVBfOAWN43g4BBZrbIzBab2WUpK116xXJuHgRmmlkZ8Dfg/6SmaL6QkrUtkzbRsHifmc0ECoEp6S6LF5hZDvAIcGOai+JVLQk0NU4lUJt/y8xGOOf2p7VU3vAl4FfOuYfN7Dzgt2ZW4JyrSnfBsoVfa2axrJUW2sbMWhKo+u9LSenSK6Z15MzsM8C/Alc6506kqGzp1ti56QQUAAvMbCuB9v0Xs2QQSCzvmzLgRefcKefcR8CHBMIt08Vybv4ZeBbAOfcu0JbARLsS4zWpufwaZrGslfYi8NXg19cC812wNzLDNXpuzGwM8ASBIMuWfg9o5Nw45w4453Kdc/2cc/0I9Cde6Zxbnp7iplQsf1P/S6BWhpnlEmh23JLKQqZJLOdmG3AxgJkNJRBm5SktpXe9CHwlOKpxInDAObcr0QfxZTOjq2etNDN7CFjunHsReJpAVX8zgc7J69NX4tSJ8dz8EOgIPBccE7PNOXdl2gqdIjGem6wU47l5FZhmZhuASmC2cy7jWztiPDffBJ40s28QGAxyY5Z8eMbMfk/gQ05usM/wAaAVgHNuDoE+xMuBzcBR4KaklCNLzreIiGQwvzYzioiIhCjMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER8T2EmIiK+pzAT8QgzGx9c76mtmXUIrhlWkO5yifiBbpoW8RAz+y6BqZDaAWXOuf9Kc5FEfEFhJuIhwbn/lhFYS63YOVeZ5iKJ+IKaGUW8pSuBeTM7EaihiUgMVDMT8RAze5HASsb9gZ7OuTvTXCQRX/DlrPkimcjMvgKccs7NM7MWQImZXeScm5/usol4nWpmIiLie+ozExER31OYiYiI7ynMRETE9xRmIiLiewozERHxPYWZiIj4nsJMRER87/8D0KpuMlQfFdsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(7, 7))\n",
    "ax = fig.add_subplot(111, xlabel='x', ylabel='y', title='Generated data and underlying model')\n",
    "ax.plot(x_out, y_out, 'x', label='sampled data')\n",
    "ax.plot(x, true_regression_line, label='true regression line', lw=2.)\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Robust Regression\n",
    "\n",
    "\n",
    "Lets see what happens if we estimate our Bayesian linear regression model using the `glm()` function as before. This function takes a [`Patsy`](http://patsy.readthedocs.org/en/latest/quickstart.html) string to describe the linear model and adds a Normal likelihood by default. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sd, x, Intercept]\n",
      "Sampling 2 chains: 100%|██████████| 5000/5000 [00:03<00:00, 1294.12draws/s]\n",
      "The acceptance probability does not match the target. It is 0.8856408031777261, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    pm.glm.GLM.from_formula('y ~ x', data)\n",
    "    trace = pm.sample(2000, cores=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate the fit, I am plotting the posterior predictive regression lines by taking regression parameters from the posterior distribution and plotting a regression line for each (this is all done inside of `plot_posterior_predictive()`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAE/CAYAAADmL9yLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd8XNWd///XmT6jNurVRZa7ZcndxsYYMCWAv/QSCBCWJRuSsNnvj+96SbJJgJRNIV/yDZCEZMPCLuBQzIaQkGRJQsdAbGNsTHGRLcnqGkkz0oymz/n9Ic1FkiVbsiU0tj/Px8MPS6M79565su/7nnLPUVprhBBCiFRgmuwCCCGEEEkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQAyil/EqpGZNdjtFQSp2plGoY8P37Sqkzj2E/a5VSe8a1cEIcIwklMSmUUrVKqWB/CLQqpR5RSqUfx/6mK6W0UspyPOXSWqdrrQ8czz4mi9Z6gdb65aNt13+eZg5432ta6zkTWjghRklCSUym/6W1TgeWAMuAr09WQY43zCb7/UKcLCSUxKTTWjcCfwQqAZRSJUqp55RSnUqp/UqpzyW3VUqtUEptU0p199ew7u3/0av9f3v7a1+n9W9/s1LqQ6VUl1Lqf5RS0wbsSyulvqSU2gfsG/DazP6vs5RS/6WUaldK1Smlvq6UMvX/7Cal1BtKqR8rpTqAu4Z+LqXUXUqpzUqpJ5VSPUqpd5RS1QN+XquUukMptQsIKKUs/Z/9mf5jHlRKfXnA9s7+GmWXUuoDYPmQ49Uqpc7p/9qslPqaUqqm/9jblVJTlFLJ87Sz/zxdM7AZsL88m4fs9ydKqfsGnJOHlFLNSqlGpdR3lFLmo/+WhRgdCSUx6ZRSU4ALgR39Lz0BNAAlwJXAvymlzu7/2U+An2itM4EK4Kn+18/o/9vd3wT3plLqEuBrwOVAPvAa8Oshh78UWAnMH6Zo9wNZwAxgHXAj8HcDfr4SOAAUAt8d4eNdAjwN5ACbgGeVUtYBP78WuAhwAwngd8BOoBRYD/xvpdT5/dve2f+ZK4Dzgc+OcEyA2/v3fSGQCdwM9Gqtk+epuv88PTnkfU8AFyqlMqAv3ICr+8sO8AgQA2YCi4HzgFuOUA4hxkZrLX/kzyf+B6gF/IAXqAN+BjiBKUAcyBiw7feAR/q/fhW4G8gbsr/pgAYsA177I/D3A743Ab3AtP7vNXD2kP1o+i64ZiACzB/ws88DL/d/fRNQf5TPeBfw1pDjNwNrB5yDmwf8fOXQfQJfBR7u//oA8KkBP/sHoGHIOT2n/+s9wCUjlEsDMwd8f+aQ/bwO3Nj/9blATf/XhUAYcA7Y9lrgpcn+9yR/Tp4/UlMSk+lSrbVbaz1Na/1FrXWQvtpRp9a6Z8B2dfTVHAD+HpgNfKSU2qqU2nCE/U8DfqKU8iqlvEAnoAbsC+DQCO/NA6z9xx6uHEd670DGNlrrBB/XAIfbxzSgJFne/jJ/jb4woP99A7cfWLahpgA1oyjfcDbRFzYA1/FxLWkafeekeUD5fgEUHONxhDiMdK6KVNME5CilMgYE01SgEUBrvQ+4tr9v53Jgs1Iql767/6EOAd/VWj9+hOONNE2+B4jSdyH+YGg5jvLegaYkv+gvcxl9n3G4fRwCDmqtZ42wr+b+/b0/oDwjOURfM9/uUZRxqKeB/6uUKgMuA04bsM8wfbXU2DHsV4ijkpqSSCla60PAFuB7SimHUqqKvtrRYwBKqeuVUvn9tQ5v/9sSQHv/3wOfMXoQ+KpSakH/e7OUUleNshxx+vqrvquUyugfIHF7shxjsFQpdXn/6Lr/Td9F/a0Rtv0b0NM/2MDZP1ihUimVHNDwVP/nye4PjH88wnF/BXxbKTVL9anqD2+AVgafp0G01u3Ay8DD9IXkh/2vNwMv0BdYmUopk1KqQim1bjQnQojRkFASqeha+vqImoDfAHdqrf/S/7NPAe8rpfz0DXr4tNY6qLXupW+wwRv9TUurtNa/AX4APKGU6qav1nDBGMrxj0CAvr6c1+lrxvqPMX6W3wLXAF3ADcDlWuvocBv2B+EGYBFwkL7a2q/oG2wBfX1pdf0/ewF49AjHvZe+EHsB6AYeoq/PDvr6uv6z/zxdPcL7NwHn8HHTXdKNgI2+2mMXsBkoPkI5hBgTpbUs8ifERFBK3UXfgILrJ7ssQpwopKYkhBAiZUgoCSGESBnSfCeEECJlSE1JCCFEypBQEkIIkTIm5OHZvLw8PX369InYtRBCiBPQ9u3bPVrr/KNtNyGhNH36dLZt2zYRuxZCCHECUkodaVosgzTfCSGESBkSSkIIIVKGhJIQQoiUIbOEi1NSNBqloaGBUCg02UUR4qTicDgoKyvDarUefeNhSCiJU1JDQwMZGRlMnz4dpdRkF0eIk4LWmo6ODhoaGigvLz+mfUjznTglhUIhcnNzJZCEGEdKKXJzc4+rBUJCSZyyJJCEGH/H+/9KQukE8eArNWyp8Qx6bUuNhwdfOdYVr0Wqueuuu/jRj3404s+fffZZPvjggxF/LsTJQELpBFFVlsVtm3YYwbSlxsNtm3ZQVZZ1lHeKk4WEkjgVSCidIFZX5PHAdYu5bdMO7n1hD7dt2sED1y1mdUXeZBftpDeRtdTvfve7zJ49m9NPP509e/YA8O///u8sX76c6upqrrjiCnp7e9myZQvPPfccGzduZNGiRdTU1Ay7nRAnOgmlE8jqijyuXzmV+17cz/Urp0ogfUImqpa6fft2nnjiCd59913+8Ic/sHXrVgAuv/xytm7dys6dO5k3bx4PPfQQq1ev5uKLL+aee+7h3XffpaKiYtjthDjRyZDwE8iWGg+PvV3Pl8+eyWNv17OqIleC6RMwsJZ6/cqpPPZ2/bjUUl977TUuu+wyXC4XABdffDEAu3fv5utf/zperxe/38/5558/7PtHu50QJxKpKZ0gknfnD1y3mNvPm2NcJIc2K4mJ8UnWUm+66SYeeOAB3nvvPe68884Rh9eOdjshTiQSSieIXQ2+QXfnybv3XQ2+SS7ZqWFoLXU8bgbOOOMMnn32WYLBID09Pfzud78DoKenh+LiYqLRKI8//rixfUZGBj09Pcb3I20nxIlMmu9OELeuqzjstdUVedJ89wkYWEtdXZHHqorccRlosmTJEq655hqqq6spKChg+fLlAHz7299m5cqV5Ofns3LlSiOIPv3pT/O5z32O++67j82bN4+4nRAnMqW1HvedLlu2TMt6SiKVffjhh8ybN29U2z74Sg1VZVmDAmhLjYddDb5hbxaEONUN9/9LKbVda73saO+VmpIQRyG1VCE+OdKnJIQQImVIKAkhhEgZEkpCCCFShoSSEEKIlCGhJIQQImVIKAlxAjvWmcOfe+45vv/9709AiY7Nyy+/zIYNG4Cjl83r9fKzn/3M+L6pqYkrr7xywss4WW655ZZxmR3+kUce4bbbbgPgwQcf5L/+67+Oe58TQYaEC3ECe/bZZ9mwYQPz588f9XtisRgXX3yxMdfeaN9jsYztcqG1RmuNyTS2e9+jlS0ZSl/84hcBKCkpYfPmzWM6xtEcy+ediH0A/OpXvzrufQx16623jvs+x4vUlISYBLW1tcydO5fPfOYzzJs3jyuvvNJYeuKvf/0rixcvZuHChdx8882Ew2EAvvKVrzB//nyqqqr453/+52GXs6ipqeFTn/oUS5cuZe3atXz00UdA3zx5t956KytXruRf/uVfBt0119bWcvbZZ1NVVcX69eupr68f9j0DPfLII1xyySWceeaZzJo1i7vvvtvY15w5c7jxxhuprKzk0KFDvPDCC5x22mksWbKEq666Cr/fD8Cf/vQn5s6dy5IlS/jv//7vQftOlq21tZXLLruM6upqqqur2bJlC1/5yleoqalh0aJFbNy4kdraWiorKwFYtWoV77//vrGvM888k23bthEIBLj55ptZsWIFixcv5re//e1hv5OXX36ZtWvXcvHFFxsh/9hjj7FixQoWLVrE5z//eeLxOAAPPfQQs2fPZsWKFXzuc58zyjv0nI103Pfff9/Yb1VVFfv27SMQCHDRRRdRXV1NZWUlTz755KDPAPDrX/+ahQsXUllZyR133GGUPT09nX/913+lurqaVatW0draesR/fwMXlDzzzDO54447WLFiBbNnz+a1114DIB6Ps3HjRpYvX05VVRW/+MUvAGhubuaMM85g0aJFVFZWGtuPm+TdzHj+Wbp0qRYilX3wwQeTevyDBw9qQL/++utaa63/7u/+Tt9zzz06GAzqsrIyvWfPHq211jfccIP+8Y9/rD0ej549e7ZOJBJaa627urq01lp/9rOf1U8//bSx37PPPlvv3btXa631W2+9pc866yxju4suukjHYjGttdYPP/yw/tKXvqS11nrDhg36kUce0Vpr/dBDD+lLLrlk2PcM9PDDD+uioiLt8Xh0b2+vXrBggd66das+ePCgVkrpN998U2utdXt7u167dq32+/1aa62///3v67vvvtv4nHv37tWJREJfddVV+qKLLjqsbFdffbX+8Y9/rLXWOhaLaa/Xqw8ePKgXLFgw6Fwmv7/33nv1N7/5Ta211k1NTXr27Nlaa62/+tWv6kcffdQ4d7NmzTLKlPTSSy9pl8ulDxw4oLXu+zeyYcMGHYlEtNZaf+ELX9D/+Z//qRsbG/W0adN0R0eHjkQi+vTTTzfKO/ScjXTc2267TT/22GNaa63D4bDu7e3Vmzdv1rfccotRHq/Xq7XWet26dXrr1q26sbFRT5kyRbe1teloNKrPOuss/Zvf/EZrrTWgn3vuOa211hs3btTf/va3h/2dJct555136nvuucfY/+2336611vr555/X69ev11pr/Ytf/MLYTygU0kuXLtUHDhzQP/rRj/R3vvMd43fS3d192LGG+/8FbNOjyA9pvhOnvHg8ftQ7y7EqLCzEbDYfcZspU6awZs0aAK6//nruu+8+zj33XMrLy5k9ezYAn/3sZ/npT3/KbbfdhsPh4O///u/ZsGGD0f8ykN/vZ8uWLVx11VXGa8laFsBVV101bJnefPNNo6Zyww03DKoVjfQegHPPPZfc3Fygbw2o119/nUsvvZRp06axatUqAN566y0++OAD43NGIhFOO+00PvroI8rLy5k1a5bx+X/5y18edowXX3zR6Pswm81kZWXR1dU1bHkArr76as477zzuvvtunnrqKaOv6YUXXuC5554zagehUIj6+vrDpsJZsWIF5eXlQF+Ndfv27cachMFgkIKCAv72t7+xbt06cnJyjHO0d+/eYc/ZSMc97bTT+O53v0tDQwOXX345s2bNYuHChfyf//N/uOOOO9iwYQNr164dVLatW7dy5plnkp+fD8BnPvMZXn31VS699FJsNpvxb2Lp0qX8+c9/HvEcDefyyy833ltbW2uUfdeuXUbTqM/nY9++fSxfvpybb76ZaDTKpZdeyqJFi8Z0rKORUBJikiiljvj9QBaLhb/97W/89a9/ZfPmzTzwwAO8+OKLg7ZJJBK43W7efffdYfeRlpY25jIe6T0jlX/ge7TWnHvuufz6178etO1IZTxepaWl5ObmsmvXLp588kkefPBBoxzPPPMMc+bMOeL7h5b9s5/9LN/73vcGbfPss8+OaR/DHXfevHmsXLmS559/ngsvvJBf/OIXnH322bzzzjv84Q9/4Otf/zrr16/nm9/85qg+t9VqNc6/2WwmFouN6n1Jdrv9sPdqrbn//vuHXafr1Vdf5fnnn+emm27i9ttv58YbbxzT8Y5E+pTEKc9sNlNSUjKuf45WSwKor6/nzTffBGDTpk2cfvrpzJkzh9raWvbv3w/Ao48+yrp16/D7/fh8Pi688EJ+/OMfs3PnTmDwchaZmZmUl5fz9NNPA30XleR2R7J69WqeeOIJAB5//PHD7tBH8uc//5nOzk6CwSDPPvusURsaaNWqVbzxxhvG5wkEAuzdu5e5c+dSW1tLTU3fkvJDQytp/fr1/PznPwf6arQ+n++wJTyGuuaaa/jhD3+Iz+ejqqoKgPPPP5/7778f3T8B9Y4dO476+davX8/mzZtpa2sDoLOzk7q6OpYvX84rr7xCV1cXsViMZ555ZsR9jHTcAwcOMGPGDL785S9zySWXsGvXLpqamnC5XFx//fVs3LiRd955Z9C+VqxYwSuvvILH4yEej/PrX/+adevWHfVzHKvzzz+fn//850SjUQD27t1LIBCgrq6OwsJCPve5z3HLLbccVs7jJaEkxCSZM2cOP/3pT5k3bx5dXV184QtfwOFw8PDDD3PVVVexcOFCTCYTt956Kz09PWzYsIGqqipOP/107r33XqBvOYt77rmHxYsXU1NTw+OPP85DDz1EdXU1CxYsGLZDf6j777+fhx9+mKqqKh599FF+8pOfjKr8K1as4IorrqCqqoorrriCZcsOnwA6Pz+fRx55hGuvvZaqqiqj6c7hcPDLX/6Siy66iCVLllBQUDDsMX7yk5/w0ksvsXDhQpYuXcoHH3xAbm4ua9asobKyko0bNx72niuvvJInnniCq6++2njtG9/4BtFolKqqKhYsWMA3vvGNo36++fPn853vfIfzzjuPqqoqzj33XJqbmyktLeVrX/saK1asYM2aNUyfPp2srKxh9zHScZ966ikqKytZtGgRu3fv5sYbb+S9994zBj/cfffdfP3rXx+0r+LiYr7//e9z1llnUV1dzdKlS7nkkkuO+jmO1S233ML8+fNZsmQJlZWVfP7znycWi/Hyyy9TXV3N4sWLefLJJ/mnf/qncT2uLF0hTkljWbpiItTW1rJhwwZ27949aWU4Ho888gjbtm3jgQcemOyiTAq/3096ejqxWIzLLruMm2++mcsuu2yyi5UyjmfpCqkpCSHEGN11113GkOjy8nIuvfTSyS7SSUNqSuKUNNk1JSFOZlJTEkIIcVKQUBJCCJEyRhVKSqn/Tyn1vlJqt1Lq10opx0QXTAghxKnnqKGklCoFvgws01pXAmbg0xNdsFTw4Cs1bKnxDHptS42HB1+pmaQSCSHEyW20zXcWwKmUsgAuoGniipQ6qsqyuG3TDiOYttR4uG3TDqrKhn8mQYjRGrr8wqli9erV47Kfm266yZj+ZryWdhCp4aihpLVuBH4E1APNgE9r/cJEFywVrK7I44HrFnPbph3c+8Iebtu0gweuW8zqirzJLpo4wR0plMY6RczRjMf+krNjH68tW7aMy34G+tWvfjWmpTtEahtN8102cAlQDpQAaUqp64fZ7h+UUtuUUtva29vHv6STZHVFHtevnMp9L+7n+pVTJZDEuBi6/MLQZRMGLscA8KMf/Yi77roLYMTlKQa66667uOGGG1izZg033HDDiMsQJBIJvvjFLzJ37lzOPfdcLrzwQqMGMn36dO644w6WLFnC008/PeJxn376aSorK6muruaMM84Ahl+aAfqWWIC+KZA2btxIZWUlCxcuNJZpePnllznzzDO58sorjaU9jvbYysClHUZawqG9vZ0rrriC5cuXs3z5ct54442x/9LEJ+No04gDVwEPDfj+RuBnR3rPybR0xRv72/Xib72g/+//fKQXf+sF/cb+9skukhgHqbB0xcDlF4YumzD05/fcc4++8847tdYjL08x0J133qmXLFmie3t7tdYjL0Pw9NNP6wsuuEDH43Hd3Nys3W63sRTGtGnT9A9+8ANjnyMdt7KyUjc0NGitP15SY7ilGbTWOi0tTWut9ebNm/U555yjY7GYbmlp0VOmTNFNTU36pZde0pmZmfrQoUM6Ho/rVatW6ddee+2wzzdwyY7k0g5aj7yEw7XXXmvsp66uTs+dO3ekX40YBxO9dEU9sEop5QKCwHrglHgyNtmHlGyyW1WRK014J6GXXx55du7jdeaZo384feCyCSM52vIUA1188cU4nU5g5GUIXn/9da666ipMJhNFRUWcddZZg/ZxzTXXHPW4a9as4aabbuLqq682lkAYbmmGgV5//XWuvfZazGYzhYWFrFu3jq1bt5KZmcmKFSsoKysDYNGiRdTW1nL66acf+eT1G2kJh7/85S+D+p26u7uNqYJEajlqKGmt31ZKbQbeAWLADuDwhU9OQrsafIMCKNnHtKvBJ6Ekxt3AJQ8sFguJRML4PhQKAUdfnmKk/ekRliH4wx/+MKp9HOm4Dz74IG+//TbPP/88S5cuZfv27Vx33XXDLs0wGsllFGDsyzCMtIRDIpHgrbfewuGQp1lS3ahG32mt79Raz9VaV2qtb9BaD39rdpK5dV3FYeGzuiKPW9dVTFKJxMniaMsvFBYW0tbWRkdHB+FwmN///vfAsS9PMdIyBGvWrOGZZ54hkUjQ2trKyy+/POz7j3TcmpoaVq5cybe+9S3y8/M5dOjQsEszDLR27VqefPJJ4vE47e3tvPrqq6xYseKon+NYnXfeedx///3G9xO1npM4frLInzjljaWJbbwMXH7hggsu4KKLLhr0c6vVyje/+U1WrFhBaWkpc+fONX72+OOP84UvfIHvfOc7RKNRPv3pT1NdXX3E491yyy3U1tayZMkStNbk5+fz7LPPcsUVV/DXv/6V+fPnM2XKFJYsWTLiMgwjHXfjxo3s27cPrTXr16+nurqaH/zgBzz66KNYrVaKior42te+Nmhfl112GW+++SbV1dUopfjhD39IUVHRsIM2xsN9993Hl770JaqqqojFYpxxxhnGAoAitciErOKUJBOyfizZt9LR0cGKFSt44403KCoqmuxiiRPY8UzIKjUlIU5xGzZswOv1EolE+MY3viGBJCaVhJIQp7iR+pGEmAwyS7gQQoiUIaEkTlkT0Z8qxKnueP9fSSiJU5LD4aCjo0OCSYhxpLWmo6PjuJ4Hkz4lcUoqKyujoaGBk2meRiFSgcPhMGbkOBYSSuKUZLVajzqljxDikyfNd0IIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZUgoCSGESBkSSkIIIVKGhJIQQoiUIaEkhBAiZYwqlJRSbqXUZqXUR0qpD5VSp010wYQQQpx6RltT+gnwJ631XKAa+HDiinTyefCVGrbUeAa9tqXGw4Ov1ExSiYQQIjUdNZSUUlnAGcBDAFrriNbaO9EFO5lUlWVx26YdRjBtqfFw26YdVJVlTXLJhBAitVhGsU050A48rJSqBrYD/6S1DkxoyU4iqyvyeOC6xdy2aQfXr5zKY2/X88B1i1ldkTfZRRNCiJQymuY7C7AE+LnWejEQAL4ydCOl1D8opbYppba1t7ePczFPfKsr8rh+5VTue3E/16+cKoEkhBDDGE0oNQANWuu3+7/fTF9IDaK1/qXWepnWell+fv54lvGksKXGw2Nv1/Pls2fy2Nv1h/UxCSGEGEUoaa1bgENKqTn9L60HPpjQUp1kkn1ID1y3mNvPm2M05UkwCSHEYKMdffePwONKqV3AIuDfJq5IJ59dDb5BfUjJPqZdDb5JLpkQQqQWpbUe950uW7ZMb9u2bdz3K4QQ4sSklNqutV52tO1kRgchhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqQMCSUhhBApQ0JJCCFEypBQEkIIkTIklIQQQqSMUYeSUsqslNqhlPr9RBZICCHEqWssNaV/Aj6cqIIIIYQQowolpVQZcBHwq4ktjhBCiFPZaGtK/w/4FyAxgWURQghxijtqKCmlNgBtWuvtR9nuH5RS25RS29rb28etgEIIIU4do6kprQEuVkrVAk8AZyulHhu6kdb6l1rrZVrrZfn5+eNcTCGEEKeCo4aS1vqrWusyrfV04NPAi1rr6ye8ZEIIIU458pySEEKIlGEZy8Za65eBlyekJEIIIU55UlMSQgiRMiSUhBAixT34Sg1bajyDXttS4+HBV2omqUQTR0JJCCFSXFVZFrdt2mEE05YaD7dt2kFVWdYkl2z8jalPSQghxCdvdUUeD1y3mNs27eD6lVN57O16HrhuMasr8ia7aONOakpCCHECWF2Rx/Urp3Lfi/u5fuXUkzKQQEJJCCFOCFtqPDz2dj1fPnsmj71df1gf03iazD4sCSUhxCnpRBo8kOxDeuC6xdx+3hyjKW+igmky+7AklIQQp6QTafDArgbfoD6kZB/TrgbfhBxvYB/WvS/sMQLxk2gyVFrrcd/psmXL9LZt28Z9v0IIMZ6SQXSyDx44Vve+sIf7XtzPl8+eye3nzTmufSmltmutlx1tO6kpCSFOWafK4IFjsaXGw3+9sZ+/X5Y74X1YA8mQcCHEKWvo4IFVFbmnZDBpramrqyO57NCuBi8//J89fG1DJVedXcX6xb2fWBOehJIQ4pQ0cPDA6oo8VlXkfqJ9J5Olq6uLgwcPkkgMXrN12rRpTJ8+HYDtvTU8vHFgH5bL6MOa6HMjfUpCiFPSg6/UUFWWNegiu6XGw64GH7euqzhhjjGz04mNAAAgAElEQVSScDjMvn37CIVCxmtKKdxuN+Xl5ZhMo+u9iUajWK3W4y7PaPuUpKYkhDglDRcKqyvyxrUmkBzhl6x9DaydjZdEIkFtbS0dHR3Ga0opbDYbs2bNwul0jmo/Xq+XAwcOEI/HUUoZr0+ZMoXCwsJxK+/RSCgJIcQEGe/pgTo6OqitrWVgC5dSimnTpjFjxoxR7SMQCLBv3z5isZjxmtYat9vNlClTqKurG7T9wO0+CRJKQghxHI7WRDdwhN+Xz545qkAKhULs3buXcDgMYNRccnJyWLx48aia3kKhEPv37x/UfAeQlpbG3LlzOXToEF6v13jd6/ViMplYsmTJqJv2JoKEkhBCHIejNdEdaYRfsumts7MT6AsfrTUOh4PZs2fjcDgOO97QEIxGozz94lbeq23jiqVlxnZ2u51Zs2YRj8fZu3evUePx+/3s3r2bV5oVp80vPyxM//zawQnv7zoSCSUhhDgOR2qiGxhQszI1Wb2HuPmeJ9l4/myqp2SjlKK8vHxUTW/xeJwDBw5g7azl7369h385fw5VZW4+aA1w71vd/OymNRQ7IzQ3NwN9Ax12796Ny+WisrLysMEKkSGjDyeiv+tYyOg7IYQYB8nZD76wpowLpikikQjPbG9gZkEaVWVu8vLymDZtGm8d7Dzi6LtEIkFdXR0ez+CHVc1mMzNmzMDtdvPKh03c+rPn+dScbP64u5WNn5pDdZmb0tJSioqKBg1UOJKjzWgRiUTo7e3F7XYf+4npN9rRdxJKQogTymQOsx4oHo9z8OBBurq62NXo454/7eGCykL+9FEnP//iRZw5v/SI79da09TUZNRsBpo2bRr5+fkAdHZ2cvDgwUGDG6xWK3+s0zy4pXHQFEDHcm6+99t3uP+5N7lmSRHXr5o+qHw2m40ZM2aQnp4+6vMyEhkSLsQpLlUu3uPtkxhmPVRbWxv19fVorY1+H5PJRHl5OW3KzU//tIP/2HgNqyvyuCpZHrvdOPdtbW0cOnTosFFzJSUlLF26FKUUiUSCAwcO4PV6qaurM0bBDTe4YUuNh6d27jisn2q4c/Olx9/huxdMY8eOHcTj8UGfa2eDl0dfrOOfLlvLkztaudRVzFRHmNbWVlpaWvD5fPh8PtauXTth53YoqSkJcZIaOmPB0O8n0/EG5kRNpNrb28v+/fuJRqODAiQ/P5+pU6cO2yw28LN0dXVx4MAB3q3vZF+bnyuXThn2/cMNy06GXHZ29qg++6DgeWw7d64vpswRYVejjx/+8SMuqCw0mvbOrq6gtLSUrq4umpqaiMfj7Kjv5HvPvM1ls2xYQ11s272HbQc8rJxVRF66DbPZhNutWbduOmeddQY5Oecd17mV5jshxDFdvD+JGtZ4BObxzGCdHDTg8/mM2g+A0+lk1qxZ2Gy2I76/p6eH/fv3H/YMT3Z2NuXl5ZjNZuDjJrqmpqZBgeZyuZg1a9aYZkqIRCLU1dXxHy++z8zCdBZN+Ti8djV4OeSHKxYVEQwGefytWp7Y2sDVS0tZXZhg3759+Hw+PB6PUbZ3D7aQk+5g3vRiysqsdGoPue5O0uxNlLg9uFztWK19w8nT0qpYvnznqMs6HAklIQQw9ov3J1XDOp7azmjfq7Wmra2NhoaGQTUfk8lkDBo4kmAwyL59+4hEIoPen5GRQUVFxaBQiUQi7N2797DngsYy+CAQCHDw4MHD9hGNRunq6sLpdJKRkTHo8yX7nHw+H7FYjHf31fHS7kYqch3safSwdsFUFs8pJyPDQV5eiFBoD/H4QbKzu0lL6yA9vQuLJT60KIMoZeeMMwIoZT7qZxh5H9KnJMQp71hmwR7vWQiOdJyxPlQKI0+kes+ls8mNdxGNRgdtX1RUxJIlS44YCpFIhH379hEMBge97nA4mDNnDna7fdDrHR0d7Ny587DBB7Nnz8blco14nGSI1NfXD+rf6ejowOPxYLVaKSkpMZ5P0lrT0dGB1+vF5XJx4MABWltb6ejooKenh2g0Sjwex+12k52dzSFvkDf21vKp6gjTilq5eH0PFvtuZpSFyM0OYjKNrRKS0C6cjtm43dXE4wEslswxvf9YSE1JiJPU8dZ4xnOBtyOVb6zB97MX95KT6GKm++O79p0NXup9Mf712rOP2PQWi8U4cOAA3d3dg1632WzMnDnzsEBJJBLU1NQMmvkAIDc3l+nTpw8780EikaC5udkYVRcIBGhtbTVmZ8jKyqKwsBCLxUIikaCzs5NYLEZubi6xWIwPP/yQ5uZmAoEA3d3dmEwm4vE4iUSCQCBAV1cXSikyMuJkZfnIze3B7e4mN9dPTm4v2e7wUc/hUL3hdCLBPJqbbXxYa2H7ASdpqpjp7iLWrj2DG264Ycz7HEqa74Q4xR1P39BEr8g6msDUWtPS0kJjY+Og95rNZioqKsjMHPmuPZFIGMO1k/tSSh3xvX6/f9jBBzNmzDhs8EE0GjWeJWppaRkUciaTifz8fHJzc9FaEwqFMJlMRm0rGAyyc+dOOjs7CYVCBAIBent7aWpqoru7G4fDgdls7q8FRUlLC1BUFKGoKMLUqQmKi6O4XB4slsAYzzr4utOpaU5Dh3KJhXKpbXbxx50mqgvzCHa24A1G+ai1lzxrnCZPBzk2zVlnrOHxxx8f87GGklASQhyTT6JPaWhg+v1+Nr+0lT1NPi5fUmo0tRUVFVFaWjpi05vWmoaGBlpaWga9bjKZmD59Orm5uYcd87QZuTQ0NNDa2srOBi/7W/1csbSMtLQ0Zs2ahcViMcq0fft2GhoaDhtKbTabKS4uJjs7G4vFYtSCkrq6utixYwcNDQ0opQgEAjQ3N+PxeIzJT4PBID6fD5vNhtudzsyZLmbMMJOXFyAry4fL1Y7d3orJFBnTuU0kTPT25tDVlUl7u5OOjnTq6jQHDkTx+2N4/b20+EJkOix0B0LkpttwZ6SRmVfCh54wy2fkMXd6KaasQn7zfhdfOL+af/3iTWMqw3AklIQQx2QiR9/FYjFqamro6ekZ9PrvP+jk7BULOWNu8bDH1FrT2tpKQ0PDoPcppSgrK6OgoGDY4Bo4+GBng5d7/rSHjZ+aw0WrKnmvJciXfvpbrp3voCL38H4gp9NJZWUlDofDmJsO+prjkjW4pqYmtNbGyLj29nYjoMxms7EMRHp6OsXFxcyYUUxubi8ul4eQpYPy0iBuVytaHwKOPNhgqGjUQkeHi+ZmOw0NJhoaTDQ322huVsRifcFsNpux2WwUFhYaIZqfn8+bNR28/mEDczLiZOpuGhsbaevqxqoUFjNGX1kMRfnClfzthd+MqWzDkVASQkwarTXNzc2Dmt6O1Hw2sDY2x6345//4M6/ta+erF86juqxvhFxDyEZzzMUXzpw57DE9Hg+1tbVAX/NaW1sbfr8fl8tlhOC+dj+PbqljdUUObx4KcOd1Z3H6vCl0dXURCATweDz09vailCIcDtPe3k4wGKS9vR2fz0dHRwc+n4/e3l4AY0CC1WqlsLCQ6dOnk5WVRU6OhalTNW63D4uliUTiIPH4AaB1zOcyELDQ3GyjqclKS4uN5mYbLS12vF4LVqudnJwcSkpKcLlcRjNla2sr7e3t9Pb2EggEiEQiRrNkNKHxBqKku+yEtJW500tZvXwRc+fOpaSkhIKCAmw2G62trXR2djJ37lxOO+20MZd7KBl9J8QpZrJmcEg+szO0iaukpIRly5aN2PTm8/moqakhHo9jU4ovLVTcfM+TXHHaHLb6s3CWZDNv4SKW9Tchfqs/tKLRKNu2bTMWpIvH48YCd/n5+SilsFgsFBQUUFhYiN/vJxqNEgqFmJWfzuqKHF54v5Vz5xfQU7ubx995le7ubtrb29Fa4/F4jCY7u92Oy+XCZrPhcrkoLCxkxYoVVFVVUVBQQHd3DZmZXVitTQSDHxGL/RGLpRGzuS8EtYYhgwFH1NZhxdPmoL3NTmOjhUOHFE1NVsBNcXExTqeTUChEV1cX3d3dJBJxenp66O7u5tChQ8RiMRKJBGazmbS0NAoLC5kxYwZlZWXMnz+fmTNnsm1/M/c++yafWZBDmdvOnkYPz+9qJmRJp6GhgV27duHxeOjp6cFqtaKUYuXKleMSSqMlNSUhThIT1ReUDLvlU7PYv38/fr+fXQ1eY8aC9PR0Zs6cafTFDBUIBNi/fz+RyOC+kczMTCoqKga9794X9vD//riLT893MSffzg9+t4vleXFe+6iZC2enUZxuwWQyMXv2bHJzc2lvb6elpQW73W7Mz5YMwUQiQU9PD8FgELPZTEdHB3/b9RFvfVhHriVKc6eXKdnpFORm4XQ6cTqdOBwOpk6dyjnnnGPsv7e3h5ycME1Nb9LTsxu7vZWMjE7c7m7s9rE1ucXj0N7eV+NpbLTS1uZgf4OdXftMpCkzeRl24vE40WgUu92O1WpFa000GiUcDmM2m8nLyyM/P5/s7GymT59OdnY2oVDIeJYqkUgMuhFI1vpqa2vZ+sEB0u0mMh1WzGZz3+i/niA+f5DcNItRowqHw8azUqeffjpPPfXUmD7ncKT5TohT0HiMmtNa09jYSHNzM0opoy/mKxfO44qzlvNeW3jYsAuHw+zbt++wBz9dLhczZ87EZrMRiURobGwctHR38ph/2PIuP//t61SVpPNeUw8XLiyiyRtmS62PNTNyuGz1vEHPC1ksFpRSOBwODh48SF1dHV1dXbS1tREKhYhGowQCAVwuF2lpaYS0hd3tET61bA4rKmcRcuRy/3NbOC03Qn66DYsljt3ejsPRSlqah8xML3l5fnJzw1itY7tORiKK1lZbfx+PjdZWO21tDtra7CQSJiM0eoJhmr1B8t3pdMWtLCjNoawgh9awhSlFeeTY+6Y+MplMhC1peKOKleV5xvx7oVCI5uZmWlpaiEQiWK1W43MHAgFCoRDBYJBwOIzJZCIjIwOLpS98EonEoLn8LBYL6enpuFwu8vLyKCwsxGazsW7dOq655poxff7hSCgJcYoay/NFA5vQkhfK5EShA2chGBh2j245wL+sdjMrZ/AUOTabjczMTLq6ug5rytNa093dzcGDBwkGg/j9fuOu3OfzUdfu5e3mODdfsJJls8po7DXxoz/sJBEOUV6YxUcfvM/KQnDb+kbFNbR30dMbZWpeGr29vWRlZZGZmTlo5gWXy0VRUd+0O16vl/fq20gzg8vcTXp6B7m5AQqKohSXRJlSEsHtDjHWBVcDARMtLX19PM3NVvYdNNPSbKa7y4Ld7iQtLQ2TM4OItjCzKAvACAN/TLGjOcjpswoozHTgjcCrtQHWz86hu6uLV3bVMLfASZpF09zexYFmD7lOhd3S93xWNBolEolgNpvJyMggLS0N6FtxNh6PY7FY6I5osrMyKC8tIjMzk3g8TntPkBB2KsuyjX6jYDBIIpEwfifGQIdYjBUrVrBp06axnZhhSCgJcQoaqaYUjUbZt2+f0UEPfUExXBNaUnJphuSDow+//CFPvfY+583L4/MXrjA61pP9MdFo1Ljw9fT0oJSis7OTmpoawuEwTqeT0tJSSkpKyMnJwev1EovFCAQC/P7ND7FEfdjjQSKRCF3BKNtrO8hzQKZNs7/FTwLFgqm5JGJR9rX4mV2UQUF2Bk6nE7/fbwxmsNvtaJ3A4eglO7uH/PxeCgtDFBdHKCwMkZExyk6eATo7zdQeMtPcbKHmoKLHl0EiWoTFko/Z/PG583QH+agtwJqF5ZRmZ1DX2snbexqZmp7AYdKYzWasVisOh4O6jgBE/JhifUEQj8cJRmKEIjEyHBYSmOiMmCjIyaIzZmHlnCmU5GQYzzNZLBbS0tJwOBx9Ief3EwqFSCQSRv9SMKap7/CT5wAVjxIIhfH6w2Q6LTjtNhwOB1lZWWRkZGCz2bDZbEbAWywWsrKyOOOMM7jgggvGfM6GklASp5yTdamG0UouU/CNswqY4ogaa/xs/NQclpXnM3PmTONueiCtNTU1Nbz//vv4fD7jtd7eXqNfY29bD49tbeaiFXN54cN2bqjKYFqW1RhAYLVa6erqIhwOk5mZSTgcNkLC4XAQjUZpaWmhrq6OcDhsDKMOBoOYTCZsNpsxOs9kMlHX7iM73UlBdt8ovXafn49aenCZIRiHOcVZmKJBQqFeCgs106crSktjFBaGyc/vxe3uxmYb2/M9ff09Flpb7TQ22/nwoIWQz4W3w0kwaMIb1rQGNVNzXBRl2IlEIkSjUUymvuY4k8mEy+WiKxBk76F2Mu2K7ghMz88ky2U3+nCi0ShKKWMQRUZGX7AGAgHa2tro6urCYrFgsVho74nQ0h0kL81MhqVviLvD4SAejxshBB+PbEwOyHC5XMZAhXBcUdutmVNeRkM8k8//r9UsnJpHOBw2ashDc0ApRVpaGhaLhYqKChYsWDDmf49Dyeg7ccqZjHV2JlOy6S3ZHPTcOw3ctjCd1bOKKCwsZPlyxbyFi3i3vouioky+u+kv5KggrkgnXq+XaDRKfWeAtu4IG1bNo6KigpKSEuNhz0AggMPh4M0P63h0Sx03rJ7GzNwYpgI/Dzz9FqunODGFuunu7jYucMmRblarFZut707c4XAYF9PkBdhsNhOJRHA6ncbQ61gsZjS/5ThMhPw+Dvl9ZGRkMLU4l7ySBMrVzfJ5JmZPaenv9+nGak2M6bxFItDSYqGpyUpjo5m2NiddXZk0ezKoaQ2Tk2bC0xNiZpGbLJedjAw7ljQTdQ1dlGZBiy9MTlYGJSV5WCwW/H4/Pp+PSCRCIBAgw26nrLiQWk83eY4YsaCfzlBg0MzhyWmD4vG40f8Vi8Uwm81GMEeicXpCUZx2K56gwpqTSVZGBna7nYKCApYtW8ayZcsoKSkhHo/j9Xrxer1G31B6ejp2u53e3l4e+vMOfve3vSzJj9Pw/lb2/C1g1M5CoZBRfq21MbAiGbqrVq0al+a70ZKakjipTPT0OJNh4GShA+9qs7KymDFjBoFAgEOHDhkzBkQiEXp6egYtp52ZmUl7xMKTu7v5/s3ns2HVAv6y8wC3/+p/uHaeizxrBJ/PR1tbG83NzUYNxufzcbClE7tJY1WJj0d2me3ELXYys/OYPWM6s8ryMJlMWCwW9jS0c7CxlbI0aG5upqenx2g2NJvNxoJ2VqsVi8WC0+mksLCQwsJCrNYoicRB0tI6yMnxk5nZRU5ON3l50TH39/T29g2pbm214/GkEY0Wk56+gFgsn/r6ho8HEITDBIPBvmmN/DEaAibmlGRTmtY3eWpTq4emzm6yHWacdiuRWIJ2f4hspwW7xYzdbictLY1YLIbX68UfDNETCGEzKyIJTYbTjq2/HyjZRJacMy85ajC5VPqcOXOoqKgg4Czg315qYuOaPNLDHl5/bx+bXtzJqlIraSpm9AElgzwWixnNglprY/RcJBLB2xumtqOXbLvC4/NTlOXAneY0AtDhcBgDHODjWlJWVhb5+fmsXr2ac84553j/GUvznTg5HEuT3ERPJDpRtNbU19fT3t5uvBYKhaivrzeGAw/V29tLd3c3OTk5pKenY7VaCYfDpKenk5uby+7du3n//fcJhUJ9zXSNbfztg1qybQk8vVFmFGSQk/bxhSx5gU3WdKZOncrUqVNxOp19Nav6erq7u4nH431LJjS28v6BQ+Q5FHaLiThm2v1RinNcZDrtZGZmGk2A0Wi0f/LRKOnpMXJyesjLC1BYGKa0NEZpaZzc3LHVegD8fhs+n5vOznRaWmx0dmbi8bhoa4sTj/eFaLImEggEjOazWCxmXLjNZjNxrfAGwzjNZoKxOG6XnXSXg964CafdgkXHCQQCxONxYgmIJTQuW9/vRCmF1WrFlZFNi07jnOWVVFWUUd/h57ltB1kzPY0M88dD1GOxmNHElnx/smkvEolQ195DpstGaX426enp2Gw2OoMxOoNxZubYCQQCdHZ20tvba9R4kkwmE1arFZPJRDCu2NMWoGpKLrOnl+LIn8qf6hP843mVLCh14/V66e7upre3d9Cw8oFDyletWsV55x3fAn/9n1Ga78SJb6xNcseyVMNk6OzsNB4AbW9vJ5HouxhnZ2djtVrx+/309vZiNps5EHJRNWsquSpoXEQOtHXT6guydFo2nZ2dvPfee/T09NDT02NcoLTWOJ1OY/ABwLTCHOLOXLbta2T+FAezC/oGK6Snp1NUVITJZDJmp+7p6WHHjh389a9/paenx7h4m83mQSO08l1WPIEYWU4z3eEIBRkOTPEoXm+AeLyR4uIIRUUJikpiTJuqKSuNkZ4+tvBJJMDrddDe7qSlxcbBehMHDkJnm5WeHvqbw4LEYj1o7TEuzANrlsmLbXL6neSoNZPJRCAcobnTT7bLgc1iwhqJ4e0Nk0jEsVstROMKR0YGCxYsMGajSDbDJQM9Ho/T6A0x02Gmu7mWd3s8OBwOVpU58MetLF9QDkB3dzdtbW14vV5juDZgnNN4PE66KULQ18Oezjajz81ut2O32+nQbrKysqiqqiIvLw+3243T6TT6toLBID09PWit2Vrfw1n5GcwrzcFi6XvGq6Coh/dqGqjIdVBWVobFYsHn89HT00MgECCRSOD3++nu7iYSiRz2fNlEk5qSSHmjbZJLheW/E4kEHo+HQ4cOUV9fz969e6mrqyMWixl3n1prXC4XZWVl5OfnG30u0DfM1+/3G30uHo+HnfvreO29WuYXOsi0mQnEFO819jCv0EGm3Ux6ejoOh2PQDAZNTU3GXXTyghaLxahv97KjrouSdDP1LR6m5zjRkYBxp5wMx2TfglLK6BtKS0vDZrMZd/PJZiOTKUFGbpS8ggizyxNMK4tSXBylpCRO/yw8oxaLKTweO01NVg4dUjQ0mGlpsdPYqAgE+oLVYrEYTVVKqcMeGk0+fJpc0jyRSBj9N8nZHgaOgnM6nfjCkJudSX5W32g2m81Gmy9AIAqLZ5aSm5uL2Wymvb3dmIrIYrFgt9v7+n8iEaM2mgzCgcGdDMnkvxGns2+4+MABCclaalZWFuXl5ZSVlRnTGEHfUPjktEfJWRuSNxFdXV10dnaSSCRIS0vD7XYb++3t7cXv9+P3+weda6UUSikyMzPJz88nLy/PmK/P6XQazYzFxcWUlZWN/T/DENJ8J04qo2mSG+/Rd8kRaJ2dnXg8Htra2ujo6DCaSzo7O/H7/UbTS3L6F601drvdWHNn3rx5OJ1OIpEInZ2dRqdyW1ub8aBjfX09TU1NmEwm0tPTMZlMxkVWKUWHP8zORh/T8jKoa/MxryiNNHNfcKSlpZFIJOjt7e1rQuqfFsfn89HS0kI4HO4bZt3dg88fxGWzYDUrYhp6w7G+GaLTnCiljGBK1iaS14e+i3mU3NxeSkqiTJmSoKwswdSpmsLCOCNM5jCiUEhR12ChtdlCc5OF1lZH/8CDBMFg1OgvSQ5tHvgMVfIPYARMMkztdruxOmtOTg7Z2dlkZWWRl5dHaWkp4XAYr9drXKi9Xi+BQMCoZcTjcSKRSN/UR/0j2ZL7j8fjWK1Wox8oeX6S4ZRsfksG+sAh28nzOWPGDBYtWkRhYSEulwuPx8NHH33EBx98MGjk40DJz59sBk2GUCQSQSlFTk4OOTk5ZGZmDqodJgMyLS2NtLQ040YnWZNKNiWGw2Gjj89ut2OxWHA4HMbXc+fOZfny5WP+/zOUhJI4aQysKf37awe5/bxZfG5txaCfHyl4EokEPp+Prq4uvF7voOac5Kiljo4OoykG+ppYOjs7jW2UUrjdbtxuN6FQiM7OTtxuNzabje7ubgKBAPn5+ZSUlJCdnY3T2XeR93g8/OUvf+HgwYNGE1jySfxkE53WmqysLNxut3HBHThMOiMjg4KCAg702nn9QBeLcuJMT9e0t7dTV1dHIBAwPmdPT49xkUnenSc/VxwzdqsFm9VsXMxD0TjhaBy7KTFgWHOM4uIwZWUJpkxJMGWKprQ0Rn7+2Pt7vF44dEhRX6+Mvw8cTNDWqlEmM3ab9bAmtWRNKD8/n/LyckpLS+nt7SUYDBo/C4VCxkOiJpOJoqIi3G434XDYqPUlf0/JkYZmsxmXy2Vc3C0WC5mZmdhsNqNpy+FwkJ2dTW5urjFgJNn3ZDabyc3NJTc3F4vFYvTrJIM8IyMDl8tlBJvP56O7u5toNHrY1D/J0XfJsHG73UataOi2A9eCStakXC6XEdyAEUTJBQEHzuQQDoexWCxGjddms5Genk52dl9/ldPpHPQ7i0ajg0Jr3rx5LF26dMy/+6EklMRJYWAT3PKpWdz3xx3c+7t3+LuVxVy2uIydDV5++KeP2Hj+HGM26WSTTVdXFx0dHWitycjIICMjA6UULS0t1NfXG4uvZWRkGHeZyafakxeCQCBgzIk28LmajIwM/H4/tbW11NXV0dvbi9frNe5gk53ZyWMnh+gmgyb50Ors2bPZUtvN1FwXZRkW2tvbeWHrB7Q31eP3duFUfc15vkCAQG8UM3Hi0RgmMygwLqaQrM183KSUvFgmBy0MrAVE4wmKSyzMmKYpLY0bAw3KyhK43WO/JrS20r98gpmGRiuNrU5aW11AOul2q3ER7omCJ2JhZlkhjUEza2cXkGFJGMH6cZOgadDAhKysLKPfJHmjkLybT4bwwJplssZYWFhIRUUFxcXFZGRk0NnZSXNzM52dncaw946ODmMZ9GTTX3IfybB0uVw4nU6jWc7v9xu1WrfbTVpamjGgIjmbQlpaGrm5uYOmRgqFQkZomkwmIwAzMjKIx+OEw2GjXycpeYOR5HK5SE9PJy0tDbvdbnzm5L/Z5E1IsnaXbOYcWOsdeN0f+Dr0jRDMzs7G7XaTnp5+xGXkx0JCSaSsYDBIZ2cnnZ2dhMPhQX0tQ+8Qn3mnkVmF6VSXuY0728ffPMCvXvyQyxYV8fsdddyyooDSNGU8r5Ns88/MzCQzMxOTyWSERigUMu4Sk00hgUCApqYmY1JtGekAACAASURBVJhtsgmot7eXjo4OWlpaCIVCxl1nIBAwHgC1Wq1kZmbicrlwu93k5OTgcrmMO2673W7crXu9XuPuM7naaCwWIxSN09bRhVnH0PEYCa3R/Rclk8mERqE1mM0mbFYLNoeTnqgm06owKYwmrqH9QX21ISgoiDNtmmbaNMW0aVBaGqe4OIrTObb/+7EYtLZajOUTksOt29sdWK2Zxt1+8vgxi5Nth3qYl6OwxEK0eEPUdfZQlGahMMtJMKY50Bmi8v9v78yj27zOO/1crARAggAIAiBBUhQpUZsl2bJd21K8JW4WZ3GSpo2bJnXTnrTxJE3bzPR0memZZjptp02baeckM2lymqRJs7htlrq1YzsnTmTXuyzZ2iWLIinuCwhiIXbgmz/Ae/lREiVS4ird5xwdccFyefF9973v9rttjUSDPmWMzJ6LXOjlwit7eaTIqFyApTFzOp04HA5lPOQCL+dDNqtKfbdIJMKWLVsIBqvNpOl0mlKppA7uM1e3SW9N5gRlUYHsM5IySzL8ar6eK5WKen8ZTpMGBVA9XFIuSEojSS/PjPl+MWP2qKRX7/V6L1q1uRpoo6RZVgzDIJVKMTk5STwen7OzM9+M5utL7tDcbreK98tdZKlUYnR0lL6+PmKxmCo7louRXHRlHue5M2O8cGqE3U0ubmz1qVCG1WpVPRcyIW82NplMhrGxMU4NjGOnjLWcV0nrbLFMOlvAYzPIZrPq6AK5gMjkdDQapa2tjWg0CkB/f786ejubzZLNZtU5NvK95S5Yliab5wsAiwUMsNpslCsCZ42NQlngsFrIF8o4bAY2gXpepQJYLdS6amaUAaxs2GBhwwZobi4RCuVobMwQDC5eTDSfh6EhG4ODNiYn68jlQmQyjUxOuhDCrsJDctGWfyOg8izSa8mVDMYKdtoiQQZSFXZuCNIVrlOGeTiWZGo6S9TnUh6e9FDMR1DIHI4M70kPxufz4ff7VRWazWZT15zP5yMajZLNZnn11Vc5cuQIw8PDqvTZfH3KMJosq7fb7XO8W3ldeTweFfKSRqRupqHVfM2bPWfZX3b+hst8T5h/brFY8Hq9yrDIAon1jjZKmgUhG/4mJydJJpMX/P58L8Z8M3m9XgKBAPX19Zc8tuD48eN0d3er3axcsNPp9BzvQwp1yp2hDIWYc0B1dXUULE4O9KfZ3tbIybEMd28O4HfM3uDJZJJ4PM7ExITKO8gFTYZhxqdSdI+luXFzCx3RMGPJHAf7k+zbEmFHRwutra1Eo1EikQi9vb0888wznDhxgr6+PhKJhDJaZu8EUKXX5pCKeQcrd8uya18u4qVSiVQ2T65QxmkXeN0u8oaFdLaE1+Oksb5a+dbYWKOaSsPhPJFInqamAo2NlUU3l6bTFvoHbPT0O5hKeClO1zMwYCWZrMHrrZYZT01NkUgkVO7CfJqpXLhDoZDKS5wcTdPU4GVzc1B5LP/2ag+n+kbY3d7A7ZurJ8vKPEpDQwOtra34/X51PcZiMSYnJ5WRSqfTKtwmDbyseHvh1BD1Dgh7nbMH3KVyTKTybG/yKhXx9vZ2brjhBjZv3kxDQ4OqXJuPfD6vcpCymEV+zt99dYDNkTp2t/jUe77WH+fMaJoP3NIKoDxnv98/p4LuekYbpeuMfD6vQmJyZwYX91TMv7Narfj9ftV8ad79yv9lbkSGTIaHh+nr61Pd+uVyWYVKzj/l0mq1qsRyNpulWCyqJLusapLGwufz0dHRQSAQAKpegTkenk6neeHIaZ545Qw/0+7FYykxlshy8Owou9r8hOo9c0p/fT4fwWAQr9dLJpOht7eXs2fPks1mq39rxcrJc0M4itOkM1mcoqyMiNztm5kNiVnn6J2Z59hut+NyudSZNObXka8r50W+BkCpbJDOFaoqAOUy4UZBNFqmo92gpaVC+wZoazPw+xd/vyYSdmKxOqanG0gkfIyOOhkedvBG3zSnB2P4XTbi+QpdTX4CtS51rQgh8Pv9RKNRdeSBnNdQKER7eztdXV1YrVaOHz9Of38/B0738+3nu7l/e4CQS/DymRFe6pliW2sDvRNp7unw0hLwzAlxyQ2JrCaT2m6VSkVdF+b8hsPhwOfz0dzczEillt979Axf+KU9c9oA/uYXdrIjaDcdiDf7WZrvg/PvD/k5OZ1OlVeR+SLJlbYeXO/ajEtmlIQQrcDXgTBgAF8yDONvL/UcbZSuDLnwypCYXNjNkiGyokbG2QEVmpJqv3Ihl695fthAVhXJBsl4PM74+DjxeFyVysrzV+RCJJ8LqPCKXIBlctVisczpvZC7apfLpcpNZXmsTNbKkNrExASxWIxUKqWS8vl8nkQiwdDQEIlEgnQ6Tf94Crso4ZwJk5dKJQplyBaKeOyzSV7pgZm9LYvFohZC1btSqf5d9poa3A7bHAMMKC9Lfi3n29xBb66Yk8j5kfMgT0iV4aV8BdxOO+1tVkKRAvXBIls3WWhtKRMJZ6hxlhZ17VQqEIs5iMVqmZrykUjUMzHh5tw5QT4/qzrgcDhoaGig7KzjmcEyn3zgTdzcGeHg2VG+/Ew3v3X/bvZubSWXyzE1NcXg4CAjIyPqCPCLbXbka8qNRu94kn97bYDOiJ/jEwXetaed+3ZEq6KuL/bzkdtb2RKpV59JXV3dnIpFeeSExWJReRVzQ/D5Yzg8mOAvnzjJO3aE+cFrw3zotlY+cMsGZVSOjhc4OpRa0oX/SuSs1kIf3WqylEapCWgyDOOgEKIOeBV4r2EYx+d7zrVulGRJr9mbkDs9YE5HtFRKhlkpEbl4mauF5Ot6PB6VLDUfYuZyuVToyemsyoyMjo4yODhIIpEgl8upqqDx8XFisZgaj/RwZNe2YRgq/CLl6WUs3el0qvJaOa5cLkcsFlOnYcqwndVqVVVLMq4vT8pMp9OcPn2anp4e1b0uE9nSuJbLZVVdJY2EzDXInwFzvCo5h+bkvuzpkcbCvHCZk7xmL0UIQakM+UoZh8VKoVzCabViERd6STK3YbPZyBaKOB1V+RmbrWrE8qUyxbKBz1OjFlbz51tb6+SGG+ppbTWora0eHufzpwg2ZHE4FltsYCGV8pHPR1TBQTbbSD7fCFTDZfKU0mAwSGdn55zwqvRMnjzUQ9Al2NTkV4fvvTE0yUgyx63tAWpqatTGQRYKyMrB89cM+X7Sg/H5fDx6NMYjL/dyR6uLB28KKUHWkYKdoWn4wM2tlwwNy/CwLJW/XF5F9rG976Zm9p+eWJGF/0rkrK5FbcaFslCjdNmWN8MwhoHhma9TQogTQBSY1yitBnJ3bDYUmUxmzsJ8Ma9hvu9lnF/2GkjxRnnjSMkPWZbqcrnUjW+329m4cSOBQEDJmJjfI5fLMTIywsDAAMPDw0pmROZDhoeHlSy9NCay4U/+LQ6HQ3klsoxVils2NDRw4403UqlUVA+Nx+OhpaWFSCTCE6emVEVbsVhkfHycg93DDCVL3LzJRywW49SpUwwPDyvjJpO90qBIY2m1WpmenlZ5FmmcZVhOeicyhi/7Pswy/ubkvzmZLY2P7OWRf/v5xtz8vexEt9vtKnRYqVRUebRcQEvlCvlyBafdit0isIhqz06Nw4bTYVMGW24MZIhSWO0k8iUcpXK158XjIzFdYVPYS9GYZkMky4aWPOFwjkikQLAxS0Ng8YfH5fN24vE6YrFa4vE6JifrSCR8QJja2up5Nzt37mTv3nasVivFYpF4PE48HlcCmzJJL6+r3t5e1af1zls3Y7fbSSaTJBIJjowV6Wpt5a6wV83zkb4xuk+c5s6uMPX19eoalfeBvF9qa2vZuHGjutZfH5jiqVMT/Mqbd/H46SThXfu4d/vVKwLMx/nSUg/f07HsC/+Vylnt7Qzy4dvalDG7XgzSYlhUH7YQoh24CXhpOQZj5sSJE6p3wWwsLmVcpHGQ5Z5ut7u6u50pQY7H4xeUV14MuQBKKZhAIKD6VORNLGPV8oAziVz4pOyHTNzLMcpdJ6BCbrICTeYsZIWPlHSRyVm3262Oe5YhDavVqsJ2UmVgdHSU7u5u1bktjbP0JJLZPL3j07QFXNTW2MmUDHrGkjR77bzgsKleF9mNXqlUlBCleb7Nf5ts6JNHMEvPTRoWc7IcZnW+zE2TQghVmCDfQ3aXy2SxDBW5XC4GJ9PYjGooTz4nNpUknZ7GKrIqt2RuRpRflwyosQqssiTX6aDGZcUQVhr9XhUCleflyFLl2loPdk+FnDXGDZsq1HnH2d5RJtzYg9s9e4DeQkmma+gd9VIshKm3t1EuR4E22tpu4m1vuwOLxcLIyMgcjxuq/TxTU1OcOXNGlSzLsmTZLDwxMaGKPYQQhMNh/H6/0jqTZfMtLS3kXZnq8RR3tNEV9nJ6NMm/HE/xmQ/fz/vuuumCvMp8PN89wecfn+Qr/+WD7O0M8oD0VJw1y7IAn+8J3d7ZwCe/dYi7uxqXbeGf7z0XGsJbD9qMq8mCCx2EELXAfuBPDcP43kV+/+vArwO0tbXd3NfXt5TjVLt+2fMhwysXCyecf/OYd7yy0U6GmOSJlTLmXywWVTmnucpK9qTU1NQowyFlOHK5nArVyV4JGSeX58dIb0l6BdLbEaJ6Oue5c+cYHh4mHo+ro6JlTgJQ+R2ZaJeLreyTkAbAHOqTRk6WXpvDXmPpAv9xYpCWOisDqTJ3bm+hwWVV7w0ofTN5Ts7k5KRqLDUXM8iciVzwzQbYLM8i8y4y3Gb+bKRxCgQCdHZ2EggEVCPqxMQE8Xic6elppqam1Gd4fl6HmX4em21WjFP+k93wgUCAUCikPExzqNCcK/N43LS22gmH89TWxvD7k9TVxamtjeF0Lk6g0jAEpVKQUilKqdRMsdiE17uTc+k2/vKJId62PcRjB8/y0VtDtNZZz3tuNeleKBQYHx9XfTSpVGpOJZ/dbsfv96ucoM1mIxgMqqZRc7WZ0+lUp79eLIF/NR7GSifzL/Z+X362m8899QYfu3PjsnhKV/o36pzSElbfCSHswL8DTxqG8bnLPX4pckonT55UC5Hc7ct8i1nHScphyJ1fMplUoS9Zdiyb4WTORC5U5xsXaXjkzl2eDCl38HKRlouiDDHJ0M7Y2JhS1pXFCbIAwGxEpdSJjP9XKhXVeNfc3EwwGKSuro5kMqm6zaXYo81mU6KdsufHZrNxtD+Ox27gtRvKwKaKBnmLk12tQfUzaUxSqRQneoc5NzpBvd2g1mG5QMBSejLSwNhsNuVxwaxUiqzikyE+s2Cm2+2uzk0ig6VcwGmfndd4eppkchqbMSvwKV//YlVR0qhJbTeHw0HFYiOWq1DvtJPMlWgN1inPSQpTynFIYVJzeXggUEdrq4WGhjTQR7ncg9U6iMs1gdW6uGOzyxULQxNeplJBTo4GCLg66fR3MDVVx/R0SW1GhBBMlmx87/Akv7yvnS2RegbSFb740zN8YIeXgHX2RFMpCyM782XItrGxkdra2jnv73Q6aWpqoqGhYUFezcVYr8d+SNbywr9cBnu9VPUtZaGDAP4BmDQM47cX8uZLYZQ+85nPMDIygmEYakcsk9ZyRy5DO9KASH0nmeuRHfY1NTWqmx5mK9BkSE+WUZt7WqTHIRsF5WJpDjuZcyC1tbW0trYSDocv2P1XKhXy+bzKTQGkUimlFCA9MVkyLT0N2aAn81XmMJJcaOUcHDs3ynd/8jIdNQUoZskWK/TFMkRqLdiY1USTBqRUMUhmS9TYIF+q4HU7qXHYyRUNrKKC1SJU+CtfLFEolKlxWJSnJkN8UlVY5pLkblyOV+qGFUslcsUSNlGVxymVyxgVA2ERWE2l0TIvJgsvzPkLeWJpqVSitrYWr9eLzWbj3FSB3pybN+3cyM/ftYvW1lae7pmmK+JVvSSZTIwjva+QKfWxszlOpdKL3T6EzTaKEIvTdMvmrKQzQcqFCGeGa3nmmAUHzbw+4OGenR38R3+On93dzv6eDH/887fy5l0b8Xq9DAwM8Nprr5HL5fjxiTGavTaCjjKpVAohBBMFCzm7n4fuu3GONI00NoFA4JK9NVfLtZCEXy8L9MW41j2wpTRKbwKeBY4A8u79Q8MwHp/vOUthlB5//HGmpqbm6IYJUZUZGRoaore3l/7+/jnlojD3gCuoGiC5w5TNeIAyFDJ0Z959m9WR7Xb7nOolGUKSORK5SEqBTbO3YbFYlKGRoTgZVvJ4PIRCIUKhEK2trWzatIn29nZyuRxnzpxhYGCAc+fOcezYMXp7e0kkEnNENmUVmyxkqFQqGAiSBQOPw0Y2X8TttCgVABn+K5fLZHJ58oUiNU47LqcTYbEzlc1T57CCEKQKZRrqarAYBtl8gel8GZ/XQyQYwOPxzBETle8tiwHM+SDpRUojbRiCXKmC02knX4IGrxuXs+pRBQIBIpEIjY2NyjvL5/P4/X5aWlrU7r+hoYFwOKzmsHfayh/8+1net93Lt3/6Or/31kZaa/s5M3iAc/FT7GlP43WP4XDEF30N5vJuioUIuUoH+wejfPCud3JqooltLVvYt6lRPe5T33yF7z13lH1NgleOn+XBG4M0Osoc6xvlscND3L8zwoaGWsLhMFu3bsXj8QCoSseGhoZlNTYLYb0sbNcyV/MZrIcNxbpvnv3MZz7D4OAgMLfQQUqLNDU1EQwGVUGAPFMmmUzS3d3NyMjInLLeYrE4R79MhtWEEEqIU8rSSC0rmA0dVfMMHnVgmqyYk1VwUj5eFlv4fL450jJSJn90dJSenh76+/vJ5XJANXcjjRmgqgjNCsAyl2Nu7DSf0SK9yOlCmXypjMfpoN4zK5kvq8kMwyBZNHAIqLEJ9d65Yoli2aDe7SSdzZPIFKmxC6azeWqsBlRmRR4lZkUC6bnKn0lPT/YiBQIBWlpa6MtYOZyq5cH7buVX79lWHfPMcd5SCFIqN2/atGmOvP/AwMDM51QgkTjDePo0Q8kz3NyRo6lxGk/tJK6a7KKvNadzAx7PNtzu2X/HxkL85iN9fPDGEF//8SE+fU8L20JuTp06pcYBcC6e4fEjo9yxtZUXBvN86N7dPHTvTiKRCMFgkBd7Jq/pXfpaZ739XVdjXNZ66HXdG6WTJ0+STCZVj0xPTw99fX2qL8ZcuCAFMuVOXVaNydyOfLzMKUkF4UgkQmdnJ8FgEKfTSTKZ5DvPHCNca2FLpF7lpU4NTjA4mWZHqHomTn19PdFoVMnXDw0NcfLkSfr7+1XF4FS2hNMisFkqquE1my9SLlfwuBwqp3S+QZJJeenVWSwWVQknHy+9ENkYKg1SrlBiMlvBV+ciWYRWvxtrpaiKF2SvidvtVuXZ09PTc6ripHdTLFfAAJvDhmumIEOG7QCV72hoaCAajbJp0ya2bdvGxo0bVY+Tx+MhGAxis9lIJpP85PBZPvOt/dwWdfL82Ul++Y4NbGqsVUZM9lHJeU+l4tTWJvH5kjQ0pGeKDEapr0/gcMx6xwvBMGx4PJtNhmcrhtFKLObm68/00WgvEKkpKz2/nokUL5+dpD9R4L7dG3n3zRtxOp3s2LGD7du3Y7VaL7uzXW8L4rXIelRfuFb7n9a9UfrkJz/J4OCg6isxqy7LhLfMG0npdyn0WV9fryrf0um08pzS6TQjIyPE43ESiQTT09PK85CvV3TU8pPTcW5tstLbP8j0VIzhiQQNLiueGju5UplsvojLOiuzLyu3zI2TuUKZyekc9W47LruNoiGYmC4RrnNgpaKq8KQBkieDSk8DUKW+5+dxzGQymWqvTxlG03mCbjtGKU86kyWVzVNjAbvdprwtc5WcLE5wuVwqX+ZwOKj1NzA8LYj43IznLdx5Qwe33LCZ+vp6pSMWCoXm9G2dX6KfSCQ4fvy4qtjrHU/y41OT/NKbd3FjRzNH+sb42pOvsLe1hvZQDZFIEadzBLt9CK93Cr8/icczhcWyuHyPxeLG7d6K1bqRnpE6vv1cjqA9ylNHDT50Wwdt9Q51Bo7b7cbv9zOQKvOPhyb5vQf28JabNjNQcPGJ7xwG4KN72+e9yS+3cC1HSOz89/zi/m6sFihXUIulNnxzWU/qC+tprItl3Rulxx57jPHxcRUOkwUHsvAhnU4rBV6puWZWZZYabIZhqCol6WXIvIcM6cm8j6z+KlYgnqvgcdhI5cuE/R4avW5SuQJDsTShOgd17po5pdHmeZQeW75UJl0EX62HZL5M1F9LYEZYUxYc2O12Fa6TBkOeFOpwOFTTomxalSFIaWBkSK1YLiOoHmUAM0ceCAFWO/Uz6gOyqVIqEMtKPmlourq66E8W+cYLfTz81p3cc2MXR4dT/NWPuvn9t3exo6luzmckDd3Jkyc5d+4cgJL/F6Kqm2YYBn6/n5fPpXFXBgi5h6ivnyIQSFHrS1PriVHrmkKIxV2HhlFHodBEsdhEsdis/hUK9cTjCXpjGf71ZJLfetfN7NveztHhFJ994hR/9N4bee+dN6q8jsS8GHz1+V4A/u4jN1/1Tb7UO9jzx/LlZ7v5s8dO8ofv3MrH7rzQEGqqrAfvYz16dYth3Rulj3/845w9e1Zpp0lJF3MBAsxKnMjKMpmLkeXL8ntzj4xUEZC9QrJoQJaEA8SzRSZyFmotRTK5HN4aO8l8maDHjn1GikY2rspKMdnnZO4pOjsSp390Ep/DoM5pUWOTckFy7GAKnc14UWbtNonMgcnSdSkTJA2NrDy02+1zDg4LBoN0dXWxbds2pSlWqVSIRCIqlAjw3YODbAp51IF5qVSKZw6f4ezIFHd2+pmYmFDjlQYoHA7jcrmo9qYZBAIGfn+SYHCaurpJ7PZhHI5hLJbEoq+DUimgenuKxWbK5RZgA3V1Lfj9gQvOirHZbLS2tvJPR+LsbvUt6kaVC9e+zgY+cV7T5dXc5Esd6z9/sXz4ng7+30/PrunQzWqyXvI068W4XClLJjO0Wsjwm9R4k4fBlUqlOd6JXJwB9RizcZDl01LyRuZBXC4X4+PjjI2NqZ4VWdqcnM4xmcri97iYFm6CQS8TyQyRoIsmn0t5NLKfZnR0lFwupxL9sp9oOptnLJnFV+cmlTewW8qUC3mlAyeLLWBWsVsWDLjdblWZFwqFKJfL1NfX4/F4VMVdXV0djY2NquFXeiWhUEidsupwOIhGoyq0aBgGXq+XXbt2IYTg6NGjHD9+XBmlbTUgUkleeeUNRkdHKRQKWK1WgoUMBw50q2IOKBMKFdm506Cm5jBe7yT79k1SKJyhXL7wCIxLUW0uDVGptOL330ggcBPJpI9KpQUhPBc8PhwO09zcfMmKtYfvabjgZ3s7g/MuRud32kvMC4V87mIWiuXo4D9fquZjd3aSypa0dM1FWE/qCxe7ni51zV6rrFmjdOjQIcbHx9XiLfuBZDm0WdhUeiyNjY1zlKgnJibU8QrSkMXjcQYHB5XHIaV88vl8VUMMC5MFCy2NfqxGCTJpRseL1DosDI0myaVs1Nisqj9HPs9qtaqm3XK5TNkwKJYMnHYbuakMdgTj01YCPjcdHR00NTURDofpTwuijT42NweU/MuJ3kH6h8aI1jtUY21TU5OSC5Ld+FKiKBAIsHv3bsLhMKVSifHxcUZHR0mn06raDy7M+8gD64QQvP7663POU5LGu60tQiCQob29loaGaXK5U2Qyx8lkTmMYc5UNspcpfDNwkqeDaMMNFAoRisUoFks7P3jdyqaQ33ScuZdkjZuT47kV2SFeauHa1VI/b0jlal73akN45sWyzmXT0jXzcHggMWe+93YG+fyHbuLwQOKSc7Rcn53m8qzZ8N2DDz5IKpWisbGRUCikJHmkMKlhGPT39zM2NqbCSdKLMguCulwuFc7KZDKq3Ftq5UkDJrXc4qkMNqsFh7WqjxafzuGyWcGoYBOCRK6Ax2HFbrWoIgXphcmKvqamJs6lDLwOaPK5VSguWRRUanzsbvEyNDREKpViLJHh2GiGu3a009HcQLLi4IcnJvjofTdy+7Z2wuGwUuOWXqH0xpLJJCMjI3Oq584/H6mxsZFDhw5x8OBBstkssVhszuNdLhe7d+/mLW+5HZ8vSSZzgunpY8Tjr5PJnKRSGVp0c6nVWk++3MZLZz1sbdpJ1L+N02NB/vyJKX73bdu5tSPEhg0b8Hqr4p+rnahdaMHCYsM/yxGO0TmlleFaD6WtBus+p/TQQw9x9OhRlROS5eFyvObQnJTPkWEzeWyxx+NRCttSJkcqPshiCXOzqwybQdXATRcreFw1RIIB5Z2VbTVUHB72dDYrHbxYLKZKzWX/iuxnisfj9KcgGvLTGfGpsU/kLcTyFu7dGqI3nucbL/Syt7OBl4YK/MH7b2NntF6dcQRV4xEMBucIY8ojxZ999llee+21OcUW8rwkq9XKtm3b2LNnD16vl7q6IrHYoRmP5wSVSg92+xBW69SiPyMhglgs7VgsGxGi+r/FshGrtZGWlhbOpCz85rdfW9BivtZLWq82t7BUi5yuvtOsV9a9Udq4cSPT09MqRCfLv+XCL4+qkE2u8sgCmXwvlUpKNVk+R/YpyZyTDAdKSSFZDBEIBPD5fEriJ51OEwhU1Qymp6fnKE6XSiVCodAcIVOpEl4ulwmHw4zn4PE3Mjz0pg42NdYyNA1fOTDOr94SZHd7iEAgwCMHBvjOKwM8eGsLH75jIxaLhebmZgqFAs899xy9vb3KIzRjsViIRqOEQiGSyeSMV5bFZpvA7Z7A4RjGbh/Gbh/Cbh/CYlmskrVAiKgyPvL/YPBmNmzYPu8x6JLFLOZrtflvKQzmanuDGs1qsy4LHcy7wK1btzI0NEQykydfMvDWOVUVmrmBU2rLyZJwKQskPRbZwGm1WkkmkwzEM4TCPnZ1ddDe3o7H4+H1s4NkbV4+/p47yefzHDlyRHkoMgQoc1jSAzEbpng8TigUIhqNqmMWzGXXALcOTPHZJ07xjhoPPzw6yh//wh3cf9t2enp6+P7TlpJzfwAAD7pJREFUL/DPTx9jV9DKd753kNSZJlr9bl4+G6MlWMcdOzezd+9eisUiP3zuIG+cG+GmljqEKOF2T+FyDWK3p2hpGQPOUamcA3KLnH07Fkub8nrc7m20tu7D79+N1Vpz2c9Lcn7Ya6G5jrUq6b9UuQWZy1jL3qBGsxZYU0bJnFQ2DIPxySnG03nC9W4qFbvyePL5vDoV0+/34/f71TEO8kC5UChEe3s7xWKRRCJBJBKhq6uLUxM5/tcPDvLmznrq3PCt/Uc5MZzgrhYHX/tad1XZoWzn3OgEXUHXHJVsn8/HPffcQ2NjIxMTE6TTaZqami7oeQFUqO3gwYOMjo5SymbpyL/BV7/9HLdt9PP603089+g3KTjqeKqvzG+9/05uiPp48WQfX3u2m3e5oMnn4t+PxWgMeXE6e0hXhrBHhvjt+0q4rT1ks2eB2XLxyoJSP24slnZstg4CgZtoaNiDx7ODmpqNWCyLuxwuVQSwmMXc/NjDAwl1SJv5dVcrHHWlifKLoQ9402guz5oL38kFyn3wHzkxMs27b26j2TtX5dvlclFfX09bW5syRo2NjWzcWFVjLpVKdHd309vbq5S2BwYGGBsbq+Z5CvDT0zG6Gms4NZ4Dq5MP7usi5BL0x7P88GSch+7dwdYmH1arlXA4jM/nU2O02Ww0NzczOTnJ4cOH1ammsVhMVb2VSiVcLpeS3Tk2EOORA0O8qSvE892TfOIdu7nzhg7+5cA5PMVJNgXdBAIWjox1U+8ZJlA/xdbIJPHkMSgPLXoehQhQU9OFz7cbj2e7ktdxOqNXfKzBYs6uWUwOxfxY+fk/fE8H5cqFhm89s9bzZhrNcrKuc0qfe+oUf/KXn2Pvhlru6goBVc/jpcE8bQ0u2uqs6nTN3olpYkUrwmKh0WlgJEfJZDI0NTURL9o5MzjObR0BFdKbmpoinU7zSm+MQ0M57tuzmR3tYb71Wpx3/cwWnjg2yn9//x7uu2kTp06d4uzZs+pAwJGREcbHx5UAaTAYpLm5WR2/3tDQQEtLC8FgUFULAlU1gWdG+Z/v2UbUmeP40Bv866Hn+dBtFVr8E9jtw2SzJyiVJhc9V1Ux0VmjI4VF7fbAFc//fMyXF7m7q5HvHxpctubQa2Hx1jklzfXOujVK8mbdQw8/eu0MH7mtjc3haiPouUSRv3/+HB+5OUyL18ap4TiPHBjkHZ1uKkaFH56c4kP7uqhMDdE9OMHz51K86+ZO2hurp8C2tbXR3NzMQN7BH33/CD9Tl+RHr57k3dt89IwleOHMOLuCdraFq822kUhEHffg8XhobW2lra1tTnLf5XIRjUapr69XRqhYLHL69HGmp09TqfRwpP91IvWjBOvGKBa7qVQWV2xQrlgpiHbeiDVxY8etbInePGOEtmC1Xhg6hOUraV0pNYG1WvRwpegSY831zro0SubdY/zkS7x8up9vHorzi7u81OSrKgjDqQI/eKWPTm+Fk0Nx7t0RZduGpuphb7FpftRX5CM/924eefoAbw2liNRWZYOk8OrpoRhPvtrN3VvDbN3QRBIn3z86hcVq4+4tEV4cLfMHD+xhd4tPqYHLMvGLMTExwNmz+ymXz1Kp9FKp9GAYfTPFBqV5n3cxqmKiVW/n4GCQ7x52E8tv4M9+/t3s29R0ReerLMfOXBqM993UzP7TE0v+Hteip6TRXO+sy+o7c1L5/z7ZR7znJFuzMfY/W+HOGzqw2+2EPTZuiNbxck+M27d1srNz9rC1tgYPm8cG+du/+nP2bIqy6/bb1XHmUDVMY45m/uLt7+Znb95CU1MTB/qT7P/Gq7xrVxN//v5dswvrznY2mxbCfH6c7u6nSSQOK+NTNUAji/47bbaGC87v8Xi24XS2IoSF57sn+OsXDrEtUsdkMkH1NPrFJdmXq9rLXCX35Wd7+PRbNy9JEYD59XUnvUZz/bKmPCUzn/rUp5iamlLHUEhZn57RKZ46Nkq0uYmxop2P7mtnc6iqXv3TU2M8eTLOQ+++mx/2GXzhwzdfdiEzh1WqvU+DPH98P2+c3c+OyCTl8llKpW6s1sXpuQE4na243Vtxu7fPMUIOR+O8z1lqD8fs1fzvD85K41xJ6GipxnapUBawLGEuHT7TaFaXdRm+M/PYY48Ri8XmVIqdnczxj6/F+eyvvZ1gOMLD3zwEVI8YODaUWLDcSqVSIpfrYXr6OCMjLzIy8jKVSg822xBW64UNqpeibFiw2jeSKnVQX7edrugedYjcy335y1abSZZjUZ4tRAjyg0NDVy1Fs1QL+2ok/XWhwdKgjbvmSln3Rkke8WDm/Bvi+e4JfuMbr7IzWs+hc1N8+q2b+didszfGc2/0c3rwde5uG2F4+CVSqaNAHw7HGBbL4k4utVhcuN1byFQ6eOx4LVuie/jnw27+2wPvZN+m5kUteiuxQM6nkfbem6LsPz2+6ovxauSNdK7q6tHGXXOlrHujtFA+99QpvvTjl3lo5xj3bkpQLveQy52mWDyD0znJYltybDYfbvf2OeXVbvc2amo2IIRFvefFKsMWs+gtxQK5WI/rdx55bUlLtxcznovtolejwu5aq+pbDbRx11wJ67LQYT5KpRL9/f2Mj59SBQaVSg+TqeNsLPfwxXdUz0zIzSjrWK3Vf5fC4Wi+wPBU8z3hSzaXXkoOZzEd+0vR3X8pRYWLeWf7T48vq4zPYo55WA1ZobUqZbTe0MoUmuVkzRklwyiTy/Vy6tRTZDInKRbPUCi8gc02hBDTcx7rdV3u1SzU1GycMTxm72crNlv9osd2ucqwldZ6W2iF3UpVtK218az2e16raOOuWU7WXPhufPz7HDv2/kU9RwgnbncX2UoHY9k27tq+D7d7Gy5X17xiolfCpcJTi/VariQuP9/7f/7pMzzfHZs3JLXSyenLhchWI1muE/RLg84paa6UdZtTmp4+ySuvbLv4L0UtfYkoGyM30h65SXk/LtdGhJg/XrcSC9KVar0tZjwXWxB+4xuvAvDRve1rIr6v8w3XNtq4a66UdWuUKpUiL73UicvVidu9jdFMK59/1sKbttzBV1/K8fkP7Vn0Inct7e7Mi/5Xn+8FqiXxa+HvupbmWaPRLC3r1ihdjKWomFrqHfxq7hjlfOzrbOAT5yWaV3PXqnfRGo1mPhZqlCwrMZir4fyk6vPdE1f0OuaKoQ/f1nbVO3eZQ5LjkUZvV8viCygWg3k+ToykLvj93s7gqhmAj9/decG8ruZ4NBrN+mPNVd+ZWcqKqaWuGFoubblLoSvINBrNtc6a9pQudernYjAv5p9+6xZlTK7U65Istfd1OZZqPhbDF/d3XzBPz3dP8MX93cv2nhqN5vplXeSUrpaVOlvoWvRYdPGCRqNZCq6pQoe1yPW0WF8Pxlej0Swv10yhw1plNUJpq8VKhyk1Gs31y5oudFjLXCzst7czeE0u2FpWRqPRrBTaKGkWJZ+kK/40Gs1yosN3mkv2XF1PYUqNRrP66EIHDaCLGTQazfKiCx00i0IXM2g0mrWANkoaYOnknDQajeZq0EZJs2yKFxqNRrNYrkmjpKVxFocuZtBoNGuFa9IorZaC93pFq3trNJq1wjXZp7QaCt4ajUajuXquSU8JdDWZRqPRrEeuWaOkq8k0Go1m/XFNGiVdTabRaDTrk2vSKOlqMo1Go1mfaJkhjUaj0Sw7WmZIo9FoNOuOBRklIcTbhRCnhBBnhBC/v9yD0mg0Gs31yWWNkhDCCnwBeAewHfhFIcT25R6YRqPRaK4/FuIp/QxwxjCMs4ZhFIDvAA8s77A0Go1Gcz2yEKMUBfpN3w/M/Eyj0Wg0miVlyQodhBC/LoQ4IIQ4MD4+vlQvq9FoNJrriIUYpUGg1fR9y8zP5mAYxpcMw7jFMIxbGhsbl2p8Go1Go7mOWIhRegXYLITYKIRwAA8Cjy7vsDQajUZzPbKg5lkhxP3A3wBW4CuGYfzpZR4/DvRd5diCgNYFujh6buZHz8386LmZHz0387NUc7PBMIzLhtGWRdFhKRBCHFhI9+/1iJ6b+dFzMz96buZHz838rPTcaEUHjUaj0awZtFHSaDQazZphLRulL632ANYwem7mR8/N/Oi5mR89N/OzonOzZnNKGo1Go7n+WMuekkaj0WiuM1bdKF1OgVwI4RRCPDLz+5eEEO0rP8rVYQFz82khxHEhxGEhxI+FEBtWY5yrwUKV64UQPyeEMIQQ101l1ULmRgjxCzPXzjEhxLdWeoyrxQLuqTYhxE+EEIdm7qv7V2Ocq4EQ4itCiDEhxNF5fi+EEP9nZu4OCyH2LMtADMNYtX9U+566gQ7AAbwObD/vMf8J+OLM1w8Cj6zmmNfY3NwLuGe+fljPzQWPqwOeAV4Eblntca+VuQE2A4cA/8z3odUe9xqamy8BD898vR3oXe1xr+D83AXsAY7O8/v7gR8CArgdeGk5xrHantJCFMgfAP5h5ut/Ad4ihBArOMbV4rJzYxjGTwzDyMx8+yJVCajrgYUq1/8J8BdAbiUHt8osZG4+BnzBMIw4gGEYYys8xtViIXNjAN6Zr+uBoRUc36piGMYzwOQlHvIA8HWjyouATwjRtNTjWG2jtBAFcvUYwzBKQAJoWJHRrS6LVWf/Naq7mOuBy87NTGih1TCMx1ZyYGuAhVw3XUCXEOI5IcSLQoi3r9joVpeFzM0fAx8WQgwAjwO/uTJDWxesyIkRtqV+Qc3KI4T4MHALcPdqj2UtIISwAJ8DfmWVh7JWsVEN4d1D1bt+Rgix0zCMqVUd1drgF4GvGYbx10KIO4BvCCFuMAyjstoDu15YbU9pIQrk6jFCCBtVlzq2IqNbXRakzi6EuA/4r8B7DMPIr9DYVpvLzU0dcAPwUyFEL9X496PXSbHDQq6bAeBRwzCKhmH0AKepGqlrnYXMza8B/wRgGMYLQA1V7TfNAtekq2W1jdJCFMgfBR6a+foDwNPGTNbtGueycyOEuAn4O6oG6XrJC8Bl5sYwjIRhGEHDMNoNw2inmm97j2EYB1ZnuCvKQu6pH1D1khBCBKmG886u5CBXiYXMzTngLQBCiG1UjZI+IK7Ko8Avz1Th3Q4kDMMYXuo3WdXwnWEYJSHEJ4EnmVUgPyaE+B/AAcMwHgX+nqoLfYZqEu7B1RvxyrHAufksUAv880ztxznDMN6zaoNeIRY4N9clC5ybJ4G3CiGOA2Xgdw3DuOajDwucm/8MfFkI8TtUix5+5TrZBCOE+DbVzUpwJqf23wE7gGEYX6SaY7sfOANkgI8uyziuk/nWaDQazTpgtcN3Go1Go9EotFHSaDQazZpBGyWNRqPRrBm0UdJoNBrNmkEbJY1Go9GsGbRR0mg0Gs2aQRsljUaj0awZtFHSaDQazZrh/wOtgIz3dd4uLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 5))\n",
    "plt.plot(x_out, y_out, 'x', label='data')\n",
    "pm.plot_posterior_predictive_glm(trace, samples=100, \n",
    "                                 label='posterior predictive regression lines')\n",
    "plt.plot(x, true_regression_line, \n",
    "         label='true regression line', lw=3., c='y')\n",
    "\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the fit is quite skewed and we have a fair amount of uncertainty in our estimate as indicated by the wide range of different posterior predictive regression lines. Why is this? The reason is that the normal distribution does not have a lot of mass in the tails and consequently, an outlier will affect the fit strongly.\n",
    "\n",
    "A Frequentist would estimate a [Robust Regression](http://en.wikipedia.org/wiki/Robust_regression) and use a non-quadratic distance measure to evaluate the fit.\n",
    "\n",
    "But what's a Bayesian to do? Since the problem is the light tails of the Normal distribution we can instead assume that our data is not normally distributed but instead distributed according to the [Student T distribution](http://en.wikipedia.org/wiki/Student%27s_t-distribution) which has heavier tails as shown next (I read about this trick in [\"The Kruschke\"](http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/), aka the puppy-book; but I think [Gelman](http://www.stat.columbia.edu/~gelman/book/) was the first to formulate this).\n",
    "\n",
    "Lets look at those two distributions to get a feel for them.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd8lfX58PHPlZM9mAkzBMKSLSDgqqsq4qjYpWjdPqW2zra/oR3Wx/7a2uqvVqt9qlWqnTjb4iraOqoVlYAIMjNASFhJgJCE7HM9f9z3SU7CSc5JclbC9X69zus+5x7nXIHkXPd3i6pijDHGdCUh1gEYY4yJf5YsjDHGBGXJwhhjTFCWLIwxxgRlycIYY0xQliyMMcYEZcnCGGNMUJYsjDHGBGXJwhhjTFCJsQ4gXLKzs3XcuHGxDsMYY/qUNWvWVKhqTrDz+k2yGDduHAUFBbEOwxhj+hQR+TSU86wayhhjTFCWLIwxxgRlycIYY0xQ/abNwhjT/zQ1NVFaWkp9fX2sQ+nzUlNTyc3NJSkpqUfXRzRZiMgi4EHAAzyuqvd2ct4XgeeA+apa4O67E7gBaAFuVdWVkYzVGBN/SktLycrKYty4cYhIrMPps1SVyspKSktLyc/P79F7RKwaSkQ8wCPA+cA04HIRmRbgvCzgNuADv33TgCXAdGAR8Cv3/Ywxx5D6+nqGDh1qiaKXRIShQ4f2qoQWyTaLBUCRqpaoaiOwHFgc4LwfAj8F/H+KxcByVW1Q1e1Akft+xphjjCWK8Ojtv2Mkk8VoYJff61J3XysRmQuMUdWXu3utMfGusdnLsne3s+SxVVz/5GreKSyPdUjG9FjMekOJSALwc+DbvXiPpSJSICIF5eX2h2jiR2OzlxueWs09L23i/ZIDvLFlP1c98SGPv1MS69BMN4kI3/5229fU/fffz9133x3VGK699lqee+65qH5mR5FMFmXAGL/Xue4+nyxgBvCWiOwATgJWiMi8EK4FQFUfU9V5qjovJyfoaHVjouaelzbyTmEF2ZnJ/PzS47nt7EkA/M/Lm3lzy/4YR2e6IyUlhRdeeIGKiooeXd/c3BzmiGIjkr2hVgOTRCQf54t+CXCF76CqVgHZvtci8hbwH6paICJ1wJ9E5OfAKGAS8GEEYzUmbNbuPMgf3t9Jkkf47bULmJk7EIDkxATuW7mV7//tE14ffwZpydZnoy9ITExk6dKlPPDAA/zoRz9qd2zHjh1cf/31VFRUkJOTw29/+1vy8vK49tprSU1N5aOPPuLUU09lwIABbN++nZKSEnbu3MkDDzzA+++/z6uvvsro0aN58cUXSUpK4p577uHFF1+krq6OU045hUcffTRu2mwilixUtVlEbgZW4nSdXaaqG0XkHqBAVVd0ce1GEXkG2AQ0AzepakukYjUmnH788mYAvnra+NZEAfC108fz0vo9bN5zmKdW7eDGMybEKMK+adwdHZs2w2PHvRcGPeemm25i1qxZ/Nd//Ve7/bfccgvXXHMN11xzDcuWLePWW2/lr3/9K+B0+33vvffweDzcfffdFBcX8+abb7Jp0yZOPvlknn/+eX72s5/x+c9/npdffplLLrmEm2++mbvuuguAq666ipdeeonPfe5z4f+heyCibRaq+oqqTlbVCar6I3ffXYEShaqe6Rtj4b7+kXvdcar6aiTjNCZcCnYcoODTgwxITeQbZ01sdyzRk8B/LzoOgGXvbqeh2e5/+ooBAwZw9dVX89BDD7Xbv2rVKq64wqkwueqqq3j33Xdbj335y1/G42krPZ5//vkkJSUxc+ZMWlpaWLRoEQAzZ85kx44dALz55puceOKJzJw5kzfeeIONGzdG+CcLnY3gNiaMnnh3OwBXnzyOzJSj/7zOmJzDlBFZbNlbzYsf7+FLJ+RGO8Q+K5QSQCTdfvvtzJ07l+uuuy6k8zMyMtq9TklJASAhIYGkpKTW6qWEhASam5upr6/nG9/4BgUFBYwZM4a77747rkau29xQxoTJgdpG/rF5HwkCV508NuA5IsJ1p44D4NmCXQHPMfFpyJAhXHrppTzxxBOt+0455RSWL18OwB//+EdOO+20Hr+/LzFkZ2dTU1MT895PHVmyMCZMVqwro6lFOX1yDsMHpHZ63gUzR5KalMAH2w+w68CRKEZoeuvb3/52u15Rv/zlL/ntb3/LrFmz+P3vf8+DDz7Y4/ceNGgQX/3qV5kxYwbnnXce8+fPD0fIYSOqGusYwmLevHlqix+ZWLrkkX+zbtchfnn5HD53/Kguz719+Uf8dd1u/mPhZG7+7KQoRdj3bN68malTp8Y6jH4j0L+niKxR1XnBrrWShTFhsLeqnnW7DpGSmMDZU4cFPf/CWU4yWblxX6RDMyYsLFkYEwavb9oLwOmTc0hPDt5v5LRJ2aQne9hQVsXuQ3WRDs+YXrNkYUwYvLbJKSGcN31ESOenJnk4Y7Iz68Drm6x0YeKfJQtjeqmusYUPth9ABD47JXgVlI/v3Le32bxmJv5ZsjCmlz7YXkljs5eZowcyJCM55OtOd0sWq4orbYCeiXuWLIzppX9tc7pSnjYpO8iZ7Q0fkMpxw7Ooa2phzacHIxGaMWFjycKYXnqv2Jcsuj/zsS/BvFdUGdaYTPj86Ec/Yvr06cyaNYvZs2fzwQfOop6/+MUvOHKk++NkMjMzexzLk08+ye7du4/af9NNNzF79mymTZtGWloas2fPZvbs2WEd2GfTfRjTC1VHmti6r5pkTwKzxwzq9vUnjh/K4+9uZ/WOAxGIzvTWqlWreOmll1i7di0pKSlUVFTQ2NgIOMniyiuvJD09PWrxPPnkk8yYMYNRo9qP43nkkUcAZxbciy66iHXr1oX9s61kYUwvrN15EFWYlTuQ1KTuTzl+wtjBAKzbdYjGZm+4wzO9tGfPHrKzs1vndcrOzmbUqFE89NBD7N69m7POOouzzjoLaF9ieO6557j22msB2L59OyeffDIzZ87ke9/7Xrv3v++++5g/fz6zZs3iBz/4AeB84U+dOpWvfvWrTJ8+nYULF1JXV8dzzz1HQUEBX/nKV5g9ezZ1ddHtcm0lC2N6wVcimDduSI+uH5KRzMRhmRTtr+GT3VXMzRsczvD6l7sHBj+nR+9b1emhhQsXcs899zB58mTOOeccLrvsMs444wxuvfVWfv7zn/Pmm2+Snd11W9Vtt93G17/+da6++urWEgDAa6+9RmFhIR9++CGqysUXX8y//vUv8vLyKCws5M9//jO/+c1vuPTSS3n++ee58sorefjhh7n//vuZNy/ogOuws5KFMb1QsMNpmF6Q3/Mv+fluoimwqqi4k5mZyZo1a3jsscfIycnhsssu48knn+zWe/z73//m8ssvB5xpzH1ee+01XnvtNebMmcPcuXPZsmULhYWFAOTn5zN79mwATjjhhNYpzGPJShbG9FBDcwvrSg8BcEJez0oWAPPHDebPH+5k9Y6DLD09XNH1Q12UACLJ4/Fw5plncuaZZzJz5kyeeuqp1iomf/4r2nWcWjzQaneqyp133snXvva1dvt37NjRWu3l+/xoVzkFEtGShYgsEpGtIlIkIncEOH6jiGwQkXUi8q6ITHP3jxOROnf/OhH5dSTjNKYnNpRW0djs5bjhWQxMT+rx+/iXLLze/jGxZ3+xdevW1rt9gHXr1jF2rDP9fFZWFtXV1a3Hhg8fzubNm/F6vfzlL39p3X/qqae2m8bc57zzzmPZsmXU1NQAUFZWxv79Xa/P3vEzoyliJQsR8QCPAOcCpcBqEVmhqpv8TvuTqv7aPf9i4OfAIvdYsarOjlR8xvTWarcKat643rUz5A5OY/iAFPYdbqCkooaJw7LCEZ4Jg5qaGm655RYOHTpEYmIiEydO5LHHHgNg6dKlLFq0iFGjRvHmm29y7733ctFFF5GTk8O8efNak8CDDz7IFVdcwU9/+lMWL17c+t4LFy5k8+bNnHzyyYBT5fWHP/yh3ep6HV177bXceOONpKWlsWrVKtLS0iL407cXsSnKReRk4G5VPc99fSeAqv6kk/MvB65W1fNFZBzwkqrOCPXzbIpyE203PLmaf27Zzy8um80lc0b36r1u+tNaXl6/hx9/fiZXnJgXpgj7PpuiPLzidYry0YD/UmCl7r52ROQmESkGfgbc6ncoX0Q+EpG3RSTg8lMislRECkSkoLzc5tcx0aOqrNvltleM7X0PphPcXlAfu+9pTLyJeW8oVX1EVScA/w34OiHvAfJUdQ7wLeBPIjIgwLWPqeo8VZ2Xk9P90bPG9NTuqnoqaxsZlJ5E7uDeVwXMynW6ha4vi00jrjHBRDJZlAFj/F7nuvs6sxy4BEBVG1S10n2+BigGJkcoTmO6bUOp86U+c/TAgD1dumvaqAEkCBTuq6a+ySYV9NdfVvOMtd7+O0YyWawGJolIvogkA0uAFf4niIj/epIXAoXu/hy3gRwRGQ9MAkoiGKsx3fJJWVuyCIf05EQmDsuk2ats2Rub3i7xKDU1lcrKSksYvaSqVFZWkpra+drwwUSsN5SqNovIzcBKwAMsU9WNInIPUKCqK4CbReQcoAk4CFzjXn46cI+INAFe4EZVtRFLJm6sD3OyAJgxeiDb9tWwofRQj+aZ6o9yc3MpLS3F2iR7LzU1ldzc3B5fH9FBear6CvBKh313+T2/rZPrngeej2RsxvSUqraVLHLDlyxmjR7IC2vL2GDtFq2SkpLIz8+PdRiGOGjgNqavKTtUx4HaRganJzF6UPj6ufsSz/pSSxYm/liyMKabfKWKGWFq3PaZNnKg08i9v8YauU3csWRhTDf57vxnhbEKCiAt2cOkYVm0eJVNew6H9b2N6S1LFsZ004YING77zHDf8xNrtzBxxpKFMd2gqmza7dz1Tx8V/mQxc7Qz9tSShYk3liyM6YbymgYqaxvJSk0My8jtjqaOdJLFVhtrYeKMJQtjumHLHudLfOqIAWFt3PaZMsJNFvuqabHpyk0csWRhTDds2etUQU0ZGZlpxAemJzFqYCr1TV4+rayNyGcY0xOWLIzphs1uycJXAoiEKW5VlO+zjIkHliyM6YbNe7pZsti0An53CTx9FVQUBj8fmDLCeW9fKcaYeGBrcBsTosZmL8Xlzupnxw0PIVmseQpe9FuiZcc7cP1rkNP1BMpWsjDxyEoWxoSopKKGphZl7NB0MlKC3GdV74OV33Wen3QT5J8BdQfhlf+AIDOoTrWShYlDliyMCdGW1vaKEEoVb98LjdUw+XxY9GP48pOQOgi2vw1F/+jy0vzsDJI9CZQerONwfVMYIjem9yxZGBOizb6eUMEat+sPw8dPO8/PudvZpg+Bz9zuPF/9eJeXJ3oSmDQ8E4BtNt7CxAlLFsaEqHWMxcggyeKT56CpFsaeCsOmtO2fcxUkJEHha1DV1aKRbZ+x2ZKFiROWLIwJUesYi2DVUL5Sxdxr2u/PyIYpF4J6YcOzXb5Fa48om1DQxImIJgsRWSQiW0WkSETuCHD8RhHZICLrRORdEZnmd+xO97qtInJeJOM0Jpiquib2HW4gNSmBMUPSOz+xphx2fQCeZJhywdHHZ3zB2W595ehjfia7va0K99X0NGRjwipiycJdQ/sR4HxgGnC5fzJw/UlVZ6rqbOBnwM/da6fhrNk9HVgE/Mq3JrcxsVC03/nSnpCTiSehi2k+tv0dUKf3U0qAEsiEs8GTArs+hJr9nb5Na5vF/mpbf9rEhUiWLBYARapaoqqNwHJgsf8Jqupfxs4AfH8Vi4HlqtqgqtuBIvf9jImJov1O28GkYZldn7jt7842UKkCICUTxp8BaNu5AYwYkEpWSiKHjjRRUdPYg4iNCa9IJovRwC6/16XuvnZE5CYRKcYpWdzazWuXikiBiBTYgu4mknzVQZO6GozX0gzb/+U8n3hu5+dNWuhsS97u9BQRYaJbuijcb43cJvZi3sCtqo+o6gTgv4HvdfPax1R1nqrOy8nJiUyAxuAsdQowsauSxZ6PoeEwDBkPg8Z0fl7+6c52+7+6HKA3eZiTmHxVYMbEUiSTRRng/xeT6+7rzHLgkh5ea0xE+b6wu6yG2u6WFHzJoDPZkyFzBNTuh/KtnZ7W2m6xz0oWJvYimSxWA5NEJF9EknEarFf4nyAik/xeXgj4ZlpbASwRkRQRyQcmAR9GMFZjOlXT0EzZoTqSPQnkddUTasc7zjZYshBpX7roxCTrEWXiSMSShao2AzcDK4HNwDOqulFE7hGRi93TbhaRjSKyDvgWcI177UbgGWAT8HfgJlVtiVSsxnSl2C1VjM/JINHTyZ+M1wulBc7zvFOCv2neSc62dHWnp/hKMYVWDWXiQNBZZ0Xkf4Fl7hd4t6jqK8ArHfbd5ff8ti6u/RHwo+5+pjHhFlJ7RcU2p71iQC4MGBn8TXPnO9suksXIgalkpiRyoLaRypoGhmamdCdsY8IqlJLFZuAxEfnAHUQX/lXqjYljha3dZrvoCeX70s+dF9qbDpsGSelwcDvUVgQ8RURaE9Q2q4oyMRY0Wajq46p6KnA1MA5YLyJ/EpGzIh2cMfGgqLXbbBcli9ZkMT+0N/Ukwqi57rUFnZ422f3MIus+a2IspDYLd/T0FPdRAXwMfEtElkcwNmPiQmEoPaF8X/ihJgtoK4V02W6R1S4GY2IllDaLB4CLgDeAH6uqr1fST0Wk835/xvQDdY0t7Dp4hMQEYVx2RuCTGqph/yZnRtmRs0J/8xDaLaz7rIkXoSyruh74nqrWBjhmU3CYfq24vAZVyM/JIKmznlBlawGFETMhKS30N/eVLMrWgrcFEo6e/szXfdYG5plYC6Ua6sqOiUJE/gmgqlURicqYONE6GC+c7RU+WSNgYJ6zol4ng/NGDUwlI9lDRU0jB2ptjigTO50mCxFJFZEhQLaIDBaRIe5jHAHmaTKmP/JV/0zssidUD9orfFrbLQKPOXXmiPINzrOqKBM7XZUsvgaswWnUXus+XwP8DXg48qEZE3shNW7vWedsR8/t/gf4rtm9rtNTfJ+9zaqiTAx12mahqg8CD4rILar6yyjGZEzcCFoNVVMO1XsgOQsG53f/A0a4DeJ713d6iq/7rJUsTCx1mixE5LOq+gZQJiJf6HhcVV+IaGTGxFh9UwufVtaSIJDfWU+ovR872xEzIKEHs+eMmOls9210pjj3HP0nOclmnzVxoKveUGfgdJf9XIBjCliyMP3a9opavArjszNISexkoca9G5yt70u/u9KHOI3cVTuhshCGTT3qFN8obksWJpa6qob6gbu9LnrhGBM/QlvDwq0+GtGN8RUdjZzlJIs96wMmi9GD0khNSmB/dQNVdU0MTEvq+WcZ00NBy80icpuIDBDH4yKyVkQWRiM4Y2KpyG0j6LLbrK9k0Z3BeB0FabdISBAm5FjpwsRWKJWs17trZS8EhgJXAfdGNCpj4kBbT6hOus021EBlESQkQs6Unn/QyOCN3L7STbElCxMjoSQLcbcXAL9zpyqXLs43pl8IWg21byOgkDMVEnsxfbivZLFnfafLrE70lSzKLVmY2AglWawRkddwksVKEckCvKG8uYgsEpGtIlIkIncEOP4tEdkkIutF5J8iMtbvWIuIrHMfKzpea0wkNbV42VFRiwitVUBH8ZUEetq47TNgFKQPhfpDULUr4CnWyG1iLZRkcQNwBzBfVY8AyUDQRm93ptpHgPOBacDlIjKtw2kfAfNUdRbwHPAzv2N1qjrbfVyMMVH0aWUtzV4ld3Aaacmd9YRyk0Vv2ivAWWbVv3QRwMTWVfNsrIWJjVDWs/AC+4BpInI6MB0YFMJ7LwCKVLVEVRuB5cDiDu/9ppuAAN4HcrsTvDGR4ruDn9hZqQL8us32MllA0HaLsUMzSEwQSg/WUd9kKwyb6AtlivKfApfhrIft+y1VoPOV5h2jAf8ydSlwYhfn3wC86vc6VUQKgGbgXlX9a7BYjQmXomDtFS1NsG+T83zEjN5/YJCSRXJiAmOHplNcXktxeQ3TR9mClSa6Qpmi/BLgOFVtiFQQInIlMA9nIKDPWFUtE5HxwBsiskFViztctxRYCpCXlxep8MwxKGiyqNgGLQ0weBykhuGLe+TxzjZIj6ji8lqK9luyMNEXSptFCdCTUUBlwBi/17nuvnZE5Bzgu8DF/glJVcvcbQnwFjCn47Wq+piqzlPVeTk5OT0I0ZjAfL2OOk0WvR253dGQCZCUAYfLoLYy4CnWyG1iKZRkcQRYJyKPishDvkcI160GJolIvogkA0uAdr2aRGQO8ChOotjvt3+wiKS4z7OBU3GqwYyJOK9XKd7vLOEyMaeTMRatI7ePD8+HJiS0VWf55pvqwOaIMrEUSjXUCjp8yYdCVZtF5GZgJeABlqnqRhG5ByhQ1RXAfUAm8KyIAOx0ez5NBR4VES9OQrtXVS1ZmKjYXVVHXVML2ZkpDEzvpFAdrp5Q/kbMgl0fOIlowmePOmwlCxNLQZOFqj4lImlAnqp2a81tVX0FeKXDvrv8np/TyXXvAWEq3xvTPW3tFZ3MNKsavjEW/lp7RG0IeHh8jhPPjspamlq8nS/zakwEhDI31OeAdcDf3dezbZCc6c+CNm4f2gn1VZCeDVkjw/fBvsTTSSN3enIiowel0dSifFp5JOA5xkRKKLcmd+OMmTgEoKrrgPERjMmYmCouDzLGwr9xW8I4803OVGeeqYpCaKwNeIpvUkOrijLRFkqyaFLVqg77Qpruw5i+qK1k0UnjdiTaKwCSUiH7OEDbxnB04EtgxTZHlImyUJLFRhG5AvCIyCQR+SXwXoTjMiZmglZDhXPkdket7RaBe0RZI7eJlVCSxS04U3w0AH8GDgO3RzIoY2KlsqaBg0eayExJZPiATmaSDceCR50Z0XUjt80RZWIllN5QR3AGzX038uEYE1u+O/YJwzKRQO0RRw7A4VJISoehE8IfgK+RO8iEgsX7a/F6lYQEWy3AREenyUJEXsSZAyogmwnW9EdFwRq397jVQ8NnQEIns9H2hi9Z7N8ELc3gaf8nOig9mezMFCpqGthdVUfu4PTwx2BMAF1VQ90P/C+wHagDfuM+aoDiLq4zps8K3l4RocZtn7RBMGgsNNdDZWHAU3zjP6zdwkRTp8lCVd9W1beBU1X1MlV90X1cAZwWvRCNiZ6gySKS7RU+IVZFWbIw0RRKA3eGO/MrACKSD3QytNWYvi3mJQsIOgOtzRFlYiGUuaG+CbwlIiU4a2+PxZ0W3Jj+pKahmT1V9SR7EhgzOO3oExprnQFzCYkwrOOij2E0ouuFkKxkYWIhlN5QfxeRScAUd9eWSK5tYUysFLtfvvnZGSQGmndp30ZAIWcKJHbSrTYc/KuhVI8aJd7WfbYGVQ3ca8uYMAtpJjJVbVDVj92HJQrTLwVvr3B7QkWyvQJgwChIHwr1h6Cq9KjDw7JSyEpJpKquiYqaxsjGYozLpq00xuXrNjshlu0V4JQkuphUUERaY7SqKBMtliyMccVFTyifIGtyT/IlC5sjykRJKFOUvyAiF4qIJRbTr/naLAIOyGtpcgbKQXjXsOhMa4+orqf9KLaShYmSUBLAr4ArgEIRuVdEjgv1zUVkkYhsFZEiEbkjwPFvicgmEVkvIv8UkbF+x64RkUL3cU2on2lMTzQ2e/n0wBFE2hYZaqd8K7Q0wuB8SB0Q+YCCrG1hc0SZaAuaLFT1H6r6FWAusAP4h4i8JyLXiUgna06CiHiAR4DzgWnA5SLSsb/hR8A8VZ0FPAf8zL12CPAD4ESctTR+ICKDu/vDGROqHZW1tHiVvCHppCYFmMYjWu0VPkMnOvNPVe1y5qPqwLrPmmgLqWpJRIYC1wL/B+cL/kGc5PF6F5ctAIpUtURVG4HlwGL/E1T1TXeiQoD3gVz3+XnA66p6QFUPup+zKKSfyJgeKOqqCgqi214BzrxTw6c7zwNUReUOTiclMYF9hxs4XN8UnZjMMS2UNou/AO8A6cDnVPViVX1aVW8BOvnLAmA0sMvvdam7rzM3AK/28FpjeiX0kdvHRykiuqyK8iQI43Os3cJETygli9+o6jRV/Ymq7gEQkRQAVZ0XjiBE5EpgHnBfN69bKiIFIlJQXl4ejlDMMcp/avKjeL2RXfCoM609orpeCKnQkoWJglCSxf8E2LcqhOvKgDF+r3Pdfe2IyDk4a2Vc7DfgL6RrVfUxVZ2nqvNycnJCCMmYwLosWVQWQcNhyBoFWcOjF9SoOc62bG3AwxOtZGGiqKv1LEbgVP2kicgcnHmhAAbgVEkFsxqY5E48WAYswelV5f8Zc4BHgUWqut/v0Ergx36N2guBO0P4TGO6zetVSiq6SBa73S/r0XOjGBVOm4UnBQ4UQ90hZ/pyP5OGWyO3iZ6u5oY6D6dROxf4ud/+auA7wd5YVZtF5GacL34PsExVN4rIPUCBqq7AqXbKBJ5157fZ6baJHBCRH+IkHIB7VPXoLiHGhEHpwTrqm7wMy0phQGqADn6+O3vfnX60eJKcdouyAtj9EUw4q93hiTYwz0RRp8lCVZ8CnhKRL6rq8z15c1V9BXilw767/J6f08W1y4BlPflcY7pj6z5nrMJxI7ICnxCrkoXvM8sKnBg6JItxQzPwJAg7DxyhvqklcJdfY8Kkq2qoK1X1D8A4EflWx+Oq+vMAlxnT52xzk8Xk4QGSRUtTW+N2tEsWAKPcBBWg3SI5MYGxQ9IpqailpLyWaaOiMFjQHLO6auD2DWPNBLICPIzpF7budUsWgZLF/k3OEqdDxkNaDMaFju48WYBVRZno6aoa6lF3+3+jF44x0ddasghUDeX7kh59QhQj8jN0EiRnQfVuqN4LWSPaHZ44LJPXNu2zRm4TcV1VQz3U1YWqemv4wzEmuppavJSU1wJtM7m242uvGBWD9gqAhAQYNRt2vOMkrikXtDvcNu2HzRFlIqur3lBrohaFMTHyaWUtjS1ecgenkZES4M+h7CNnG4vGbZ/Rc51ksfvoZOFbj3vbPitZmMgK1hvKmH5t617nSzZge0VDjdNmIZ7ojtzuyFcFVrr6qEOThmeSILC9opaG5hZSEq1HlImMrqqhfqGqt4vIi4B2PK6qF0c0MmOiYGtX7RXZoBmOAAAgAElEQVS714K2wMjZkBzKONQIyV3gbEvXgLfFmWTQlZrkYdzQDEoqainaX8P0UQNjFKTp77qqhvq9u70/GoEYEwvbuuoJtesDZzvmxChGFMCAkTAoDw7thP2bYcSMdoePG5FFSUUtW/dWW7IwEdNp11lVXeNu38aZC+ogcABY5e4zps/z9YTyTZ3Rzq4Pne2YBVGMqBO+hOVLYH58gwl9XYCNiYRQpii/ECgGHgIeBopE5PxIB2ZMpNU3tbCjspYEgQkd17Hwev2SRYxLFv4x+GLyM8VNFlssWZgI6qoayud/gbNUtQhARCYAL9O29oQxfVJxeQ1edZZRPWqqjMpCqD/kzDQ7MDfwG0STr3QTsGThjNy2koWJpFCmKK/2JQpXCc5kgsb0ab4qqK7bKxaAyNHHo23YdEjKgIPboWZ/u0N5Q9JJS/Kw93A9VUds1TwTGZ0mCxH5goh8ASgQkVdE5FoRuQZ4kbbZYI3ps3zdZgPOCRUvjds+nkTIdbvQdqiK8iQIk902ly17D0c7MnOM6Kpk8Tn3kQrsA84AzgTKgbSIR2ZMhHU5gaDvCzkvTpIF+LVbvH/UodZG7n1W6DeR0dWgvOuiGYgx0dY6geCIDo3bNeVQsQ0S02D4zBhE1okxJznbT49eqNLXbrF5jyULExlBG7hFJBW4AZiOU8oAQFWvj2BcxkRUTUMzZYfqSPYkMHZoRvuDO95xtnknQmJy9IPrTN6Jzmjy3R9BQzWktJWIprR2n7VqKBMZoTRw/x4YgbNy3ts4K+eFdPsiIotEZKuIFInIHQGOny4ia0WkWUS+1OFYi4iscx8rQvk8Y0Ll+1KdMCyTJE+HP4Md7zrbcadFOaogUrKcNTW0BXa2r4ryVUNt21eD6lETLhjTa6Eki4mq+n2g1p0v6kIgaEWuiHiAR4DzgWnA5SIyrcNpO3GWbv1TgLeoU9XZ7sOmFjFhtWm3kyymjgzQXuErWcRbsgDId2PyxejKzkwhOzOZmoZmSg/WxSAw09+Fkix8ffEOicgMYCAwLITrFgBFqlqiqo3AcmCx/wmqukNV1wPebsRsTK9tcuv2p43ssLpc9T6nvSIpPbYzzXZm3Gec7fZ3jjo0xcZbmAgKJVk8JiKDge8DK4BNwE9DuG40sMvvdam7L1SpIlIgIu+LyCXduM6YoDbvcUoWRyWL1vaKk8CTFOWoQjDmJEhIhD3roL59+4T1iDKRFLSBW1Ufd5++DYyPbDjtjFXVMhEZD7whIhtUtdj/BBFZCiwFyMvLi2Jopi9r8WrreISpnSWLeKyCAkjJdBZiKv0Qdq6Cyee1HjrOpv0wERTK3FBDReSXbkP0GhH5hYgMDeG9y4Axfq9z3X0hUdUyd1sCvAXMCXDOY6o6T1Xn5eTkhPrW5hi3o7KW+iYvIwemMjijQ2+neG3c9tdJu0XrHFF7rEeUCb9QqqGWA/uBLwJfAiqAp0O4bjUwSUTyRSQZWIJTjRWUiAwWkRT3eTZwKk71lzG95mvcPqoK6tAuqCyC5ExnKdN45UtkxW+12z15eBaeBKG4vIa6xpbox2X6tVCSxUhV/aGqbncf/wMMD3aRqjYDNwMrgc3AM6q6UUTuEZGLAURkvoiUAl8GHhWRje7lU3GmGfkYeBO4V1UtWZiw8LVXHFUFVfQPZzv+zPhsr/DJO9kZMLhvA1Tvbd2dmuRh0rBMvGrTfpjwC2XW2ddEZAnwjPv6SzgJIChVfQV4pcO+u/yer8apnup43XtAHA2dNf3JJl/j9qhOksXEs6McUTclpTpVUYWvQdE/Yc5XWg9NGzWALXur+WT3YebkDY5hkKa/6WoiwWoROQx8FWccRKP7WI7bqGxMXxSwZNHSBCXuml4Tz4lBVN008Vxn60twrhnuSnmbdldFOyLTz3U1N1SA0UrG9G2VNQ3sO9xAerKHsUP81tXe9QE0VkP2cc4SpvHOV/opfgNamp1ZaYHpbmnpkzKrhjLhFUqbBSJysYjc7z4uinRQxkTKJ7vbShUJCX7rVPju0CedG4OoemDoBBgy3lmgaffa1t2+qrWte6tparGxriZ8Quk6ey9wG05vpE3AbSLyk0gHZkwkbCg9BMCs3IHtDxT2kfYKf77qssLXW3dlpSYxbmg6jS1eivbXxCgw0x+FUrK4ADhXVZep6jJgEc78UMb0ORvKnLr8maP9kkVVmdOzKCkd8k6JUWQ94Gu3KGzf32S6227xSZm1W5jwCakaChjk93xgp2cZE+c2lDpfoO1KFltedrYTz3Z6GvUV+ac7S63u+RgO7WzdPX20UxW1cbe1W5jwCSVZ/AT4SESeFJGngDXAjyIbljHhV1HTwO6qejKSPeRn+y14tNkdKzq1j01unJQKkxc6zze/1LrbV7LYaD2iTBh1mSxERIB3gZOAF4DngZNVNZQR3MbEFV8V1PRRA/H4GrdrK+HTfzuT801aGMPoemiK299k84utu2aMaitZtHhtbQsTHl0mC3VWUXlFVfeo6gr3sbera4yJV74qqJn+VVDbXgX1Qv4ZkDaokyvj2KSF4El2JhWs2Q/A0MwUxgxJ40hjC4X7bVJBEx6hVEOtFZH5EY/EmAhbH6i9wndHPvVzMYgoDFIHwPizAG1rewGOz3US38e7DsUoMNPfhJIsTgTeF5FiEVkvIhtEZH2kAzMm3Hy9g2b4ekLVH4biNwGBKX24g58v0W1um6dz9hgnWazbZe0WJjxCmRvqvOCnGBPf9lfXs/dwPZkpieQPzXB2bl4BLQ0w9jOQGcrij3FqyoXw0u1Q8pZTFZU5zC9ZWMnChEdXc0OlisjtwH/ijK0oU9VPfY+oRWhMGHy00/nSnDl6YNvI7fVuP43jL4tRVGGSPsRpu1AvbHgOaGvE37avmiONzTEO0PQHXVVDPQXMAzYA5wP/G5WIjImAtTsPAjB3rNuIXVXmrGPtSel7XWYDmXWps3UTYFqyh+OGZ9HiVRtvYcKiq2QxTVWvVNVHcaYlj+Olw4zp2kefOiWLub5puzc8Cygcd37f7AXV0eTzIWWAszZ3+VYAjh9jjdwmfLpKFk2+J+5CRsb0SU0tXtaXOV+Yc/IGg2pbFdSsPl4F5ZOUCtMWO8/dn22Omyw+smRhwqCrZHG8iBx2H9XALN9zd52LoERkkYhsFZEiEbkjwPHT3bW9m0XkSx2OXSMihe7jmu79WMa02bznMPVNXsZnZzAkIxl2fwT7N0HakL6xdkWojl/ibD9eDi3NrSWLdTstWZje6zRZqKpHVQe4jyxVTfR7PqCz63xExAM8gtPeMQ24XESmdThtJ3AtzuJK/tcOAX6A0213AfADEbFlv0yPrP3Uaa9oXTmu4AlnO/sKSEyOUVQRkHcKDJkAh8ugcCUTh2WSlZpI2aE69lTVxTo608eFOpFgTywAilS1RFV9K+wt9j9BVXeo6nqg48T75wGvq+oBVT0IvI7TI8uYblvj3lnPHTsI6g7ChuedA/Ouj2FUEZCQ0PYzrX4CT4JwwlgnQa7ecTCGgZn+IJLJYjSwy+91qbsvbNeKyFIRKRCRgvLy8h4Havo3X8libt5gp4qmuQ7Gn+ksINTfzL4CElOh+J9woIT544YAsHr7gRgHZvq6SCaLiFPVx1R1nqrOy8nJiXU4Jg7tP1xP2aE6MlMSmTwsEwqWOQfm3RDbwCIlfQhM/4LzvOC3bclihyUL0zuRTBZlwBi/17nuvkhfa0yrD90vyTl5g/BsfwMqtkHWSDjughhHFkHz3US49nfMGuYh2ZPA1n3VVNU1dX2dMV2IZLJYDUwSkXwRSQaWACuCXOOzElgoIoPdhu2F7j5juuX9kkoATho/FN79hbPzxK+BJ5SZbvqo3Hkw5kSoP0Tq+j8yK3cgqm3Vccb0RMSShTs242acL/nNwDOqulFE7hGRiwFEZL6IlAJfBh4VkY3utQeAH+IknNXAPe4+Y7plVbGTLD6btRN2vOMMXOtvDduBnHq7s131MCeOzQLaSlnG9EREb69U9RXglQ777vJ7vhqniinQtcuAZZGMz/Rv+6vrKS6vJS3Jw3FFvraK6yD1GFgZePIiyJkC5Vu4UN7lEfIosGRheqFPN3Ab05UPSpwvx4tHV5Ow5SVnkaCTvhHjqKIkIQFOuRWA4wqfwCNePt5VZZMKmh6zZGH6rVVue8XS5j8BCnOuhKwRsQ0qmmZ+GQbl4TlQyC1D19DY4rXxFqbHLFmYfuv9kkpmSTETKt5wxh6c/p+xDim6EpPhzO8AcF3jcpJp4t1CG49kesaShemX9h+up6S8lv9OftbZsWApDBgV26BiYdalkDOVgY17WOJ5g3eLKmMdkemjLFmYfundogpOS1jPqbLe6QH1mW/GOqTYSPDAZ78HwG2Jf6Fszx7KqxtiHJTpiyxZmH7p3c1l3J34lPPitG85I5uPVVMuhLyTGSqH+Wbic7xXXBHriEwfZMnC9DstXiWv6CkmJOyhcdAEOOmmWIcUWyJwwX14SeBqz2sUb/gg1hGZPsiShel3Nm7eyFe9zlrUSRfd37+mIe+pETM5OP1qPKJ8tuRnqLcl1hGZPsaShelfVBnw2rfIkAY2DjoLmfjZWEcUNwZfeDeVDGK2bmbvP34Z63BMH2PJwvQva37LuKoPOKCZHDzzx7GOJq4kpA/mlbFO9+Hs938MlcUxjsj0JZYsTP9xoATvSqfnzz3eGzhh2nExDij+jD75y/yl5VSSvA3w169Di43oNqGxZGH6h6Z6ePZaEppqeanlJA7mX0RasifWUcWdUyZk81OuZ68Ohl0fwJv/E+uQTB9hycL0DyvvhD0fsy9xJN9puoHzZxxD03p0Q2qSh+Mnj+PWxpvxkgDvPgDbbPZ/E5wlC9P3ffQHKFiGelJYWncLtQmZLJxuyaIz504bwYc6lWcHXuvseGEpVBTFNCYT/yxZmL5t+7/gxdsAKJh2Jx+3jOPk8UMZkmHdZTvz2SnDSBD4fsXZNE08H+oPwZ8uhSM2hbnpXESThYgsEpGtIlIkIncEOJ4iIk+7xz8QkXHu/nEiUici69zHryMZp+mjyrfC01eBtxlOvpmHq04F4IKZI2McWHwbkpHMqROzaWwRXsj/AYyYCQeKYflXoKku1uGZOBWxZCEiHuAR4HxgGnC5iEzrcNoNwEFVnQg8APzU71ixqs52HzdGKk7TRx0ogacudu6Kj7uAQ6d+j38XVeBJEM6bPjzW0cW9z88ZDcBzGw7C5U8765LvfA+euQaaG2McnYlHkSxZLACKVLVEVRuB5cDiDucsBtwJfHgOOFtEJIIxmf7g0C54ajHU7IWxn4EvPsErG8tp9ionjR/C0MyUWEcY986bPoK0JA+rdxxkV8tguOovkDYEClfCC/8HWppiHaKJM5FMFqOBXX6vS919Ac9x1+yuAoa6x/JF5CMReVtETotgnKYvqSiEZYugaieMngdXLIfkdJ4ucH7Vvjg34Cq9poOMlEQWuiWwv60rg2FTnYSRMgA2/Q3+vAQaa2McpYkn8drAvQfIU9U5wLeAP4nIgI4nichSESkQkYLyclvUpd8rW+MkisOlMOYkuPJ5SMliy97DfLzrEFmpiZw/w9orQnWJWxX1wkdlqCqMmg1X/w3Sh0LRP5xqPmv0Nq5IJosyYIzf61x3X8BzRCQRGAhUqmqDqlYCqOoaoBiY3PEDVPUxVZ2nqvNycnIi8COYuLH+WVh2PhypgAmfhategLRBADy92ilVLJ49ygbidcNpE7PJzkyhpLyWD7e7SWH0XLj+NRiYB2UFsOw861ZrgMgmi9XAJBHJF5FkYAmwosM5K4Br3OdfAt5QVRWRHLeBHBEZD0wCSiIYq4lX3hb4x91uPXoDnHCt0yCbnAFAQ3MLf/nIuQdZMj8vdnH2QYmeBC5f4NzP/W7Vp20HsifCDSth2DSo2Aa/OQs2vxijKE28iFiycNsgbgZWApuBZ1R1o4jcIyIXu6c9AQwVkSKc6iZf99rTgfUisg6n4ftGVbXy8LHm0C743WJnlLF44IL74aJftJty/JUNezh0pInpowYwY/TAGAbbN33lxLF4EoS/b9zLniq/brMDRsH1K2HqxdBwGJ6+El6/y3pKHcNEVWMdQ1jMmzdPCwoKYh2GCZf1z8LL34aGKsgYBl98HMaf0e4UVeX8B99hy95q7v3CTJYssJJFT9z0x7W8vGEPN581kf84r8Pki6qw6mF4/QegLTB8Jnz+/zljM0y/ICJrVHVesPPitYHbHKsOfgp/vsKpdmqoguMuhG+sOipRALy1rZwte6sZlpXC5+d27GhnQnXNKeMA+POHO6lv6rAokgiccgtc+xIMGgv7NsBjZ8Jb9zqTN5pjhiULEx+aG+Bf98EjJ8LWlyE5Ez73ECz5I2RkB7zk0bed9Riu/0w+KYnWsN1T88cNZsboAVTWNvLHD3YGPmnsKfD192D+V50R82/9BH51Imx+ySl9mH7PkoWJrZZmWPt7eGguvPE/0FwHM74INxfACdc4d7YBrN15kPdLDpCVksgVJ1r1U2+ICLef7XQ2/H9vFXGksZM1LlIy4cL74ZqXIGcKHNwBT38Ffncx7FodvYBNTFiyMLHR0gQfPw2/OglW3OyMnRg2zenn/6VlMKDz8RKqyr2vbAHgqpPHMiA1KVpR91tnTx3G8WMGUVHTyO/9e0YFkn8a3Phvp8NB2mBnMscnzoHffwF2fRidgE3UWbIw0VV3CP79EDx4PPxlKVQWwuBx8IXfwI3vwvgzg77Fyo17+XDHAYZkJPO1MyZEOuJjgojwrXOd0sWv3y6mqi7IdB+eRFjwVbhlLZz2bafasPif8MS58MRC2PCcTRnSz1iyMJGnCp++B3+5Ef53Crz+fThcBtnHwcUPO1VOsy6FhODtDo3NXn7yqlOq+OY5kxiYZqWKcDl9UjYL8odw8EgT963cEtpF6UPg7Lvg9g1O0kgZ6KzA9/wN8MAMeOunThdo0+dZ11kTORWFsPGvsH45VPqNAs4/A06+CSaeCwndu195+I1C7n9tGxNyMlh5++kkeux+J5y27q3mwofeoUWVF75+CnPyBnfvDRpqYP3T8OFjUO6XcMZ+xrkhmLa4deS9iQ+hdp21ZGHCRxX2bYQtL8Omv8L+TW3HskbC7K/AnCthSH6P3v6TsioueeTfNHuV39+wgNMm2RQvkXDvq1v49dvFTB05gBU3n0pSTxKyKmx/G9Y8CVtfhWa3m60n2blZOO585zFgVFhjN91nycJER005lLwJxW84j5p9bcdSBsKUC2H6JTDhbKeeu4fqm1pY/PC/2bqvmqtPHss9i2eEIXgTyJHGZhY+8C9KD9Zxw2fy+f5FHZeh6ab6w850IeufdhrD8fvOGTkbJpwF+ac7k0Mmp/fus0y3WbIw4afqdJfc9SHset/Z7vuk/TlZI2Hi2TDtEucOMrH3y5uqKne+sIHlq3eRn53By7d+hvTkniceE9yaTw9y2aOraPYqv75yLovCNZtvTbmzZsbWV52bi6YjbccSkiB3How7DfJOhFFznTYRE1GWLEzvqMKhnbB3g/tYD6UFULu//XmeFBh3qjMT7ISznXURwrx+1aNvF/OTV7eQkpjAszeezKxcq/OOhife3c4PX9pEZkoiT3/tJKaPCvPcW011sONdp7Sx/V+w52PalToAhkyA0Sc4j5GznO7V1uYRVpYsTGi8XqdnUmUhVBY7jdL7NjoJoqHq6PPTh8KYE2HMAmc7ag4kpUUsvL9+VMY3n1mHKvzqK3Ntfe0oUlVu+fNHvLR+D0MyknnmaycxcVhW5D6w7pDTa27Hu8706Hs+bmvr8DdgtHNTMmwaDJ8O2ZNhyHhLIj1kycK0aayFqlKnC2OV+6gscpJDZbEzajqQ9Gznbm7ETBgxC0YeD0Mnhr3k0Jnfr9rBXSs2ogr/vWgKXz/TxlREW0NzC0t/t4a3t5UzLCuF3143P/wljM60NDk3LmVroGytU+VZviVwAgHnRmbI+PaPweOcRvTMEb1qM+vPLFkcC1qancWAavZBzX73sc95VJU61UhVpVAXZHb3jGGQPQmGToChk5y7thEzIXN41BKDv8ZmL/et3MJv3tkOwB3nT+FGG3wXM3WNLVz35Ie8X3KAtCQPD1x2fPjaMLrL2+K0m+3b6PS227fRueE5UNL5TQ+AJDjtaQNGOSWTgbnONmsEZA5z/gYysp0R6TH4nY8lSxZ9TUsz1Fc5X+x1B9seRzq+roTacich1FZwVB1vIJ5k549jYK6zAtrAXOeuK3uiUyccR8X3wn3VfPvZj1lfWoUnQfjh4hk291McqG9q4TsvbOCF1oWmxnDnBVPjZ1CkKlTvdZLGgRI44Jaaq3bB4d3te+l1JSEJMnKcxJE5rO156iAnkaS529RBbc9TBnZ7vFA8sWQRaapOcbipzu9RCw3VzsCkhmporA7y2t3WHw7cPhCUuL/Uw51f6szhzi945nAYOBoGjnEeGTlx/8u8/3A9v3yjiD99uJMWrzJ6UBoPXT6HE8Z2c1CYiRhV5fF3tnPfyq00tnjJyUrhpjMnsGRBHqlJcT7rb3MjVO9x2ueqypzt4TJnX025cwNWW+4s9NRtAqkD2pJISpYz/UlKprvNCrzP/3VyBiSmQlJ61KvL4iJZiMgi4EHAAzyuqvd2OJ4C/A44AagELlPVHe6xO4EbgBbgVlVd2dVn9ThZ7NvojDJuOuL35X/Emau/dd8Rd3+HfWElbXcqaYMhbYjf88FOF0Lfc19CSM/u0/WwzS1ePth+gGcLdvHS+j00e5UEgcsX5PFfi6bEz12raadofzV3PL+Bgk8PApCTlcKXTsjli3NzmTgsM8bR9VJTvZs49jsl95r9zuv6Q04DfN1Bv+eHnOc9SjBdSEh0kkZiqtN5xPdITOvw2k0uSe72M9+ExJRuf1zMk4W7hvY24FygFGdN7stVdZPfOd8AZqnqjSKyBPi8ql4mItOAPwMLgFHAP4DJqtrS8XN8epwsPnkBnruu+9eB023U9x+VmOrcHbS7i8jq/HXHfamD4v7uv7dqGpop3FfNJ2VVfLD9AP8uquDgEWeyuQSBhdNG8K2Fk5k8PII9bkxYeL3Ka5v28dA/C9m0p+3LMj87g9MmZTN7zCCmjxrIhJyM/j8li68Kud5NJq21BjVttQeNNU5Hk4Yat4ahxu+car/aiSOg3p7F8d19zvdRN4WaLCJ5W7oAKFLVEjeg5cBiwG8OCBYDd7vPnwMeFhFx9y9X1QZgu7tG9wJgVbiDbBg6heZT/hMS0/AmOtlbE1PxJqWiiU4SUHefJrWd4/Wktk58p37tBv651z8N+yfldvlZgXrQunq/XaG/X/t9/j9ZKO8R5PPaXRf4nGavcqSxmfqmFuoavdQ1tVDX2Mzh+mbKqxtaH6UHj7C76uheLONzMrhgxkiWLBhD7mAbvdtXJCQIi2aM4Lzpwyn49CDPFuzi75/sZXtFLdsravmdO815cmICuYPSGD04jVED0xiUkcSgtGQGpiUxIC2R1EQPyYkJrY8U9yEiCJAggggI7lYIcMzdF9N2abdKKSW3d2+jCt4mpKkOmuuR5iNIUz3S7L5uqkOa6zq8rkeajpDmSSaS/wSRTBajAf/pJkuBEzs7R1WbRaQKGOruf7/DtRFZN/Mf5YO56Y05nRxtBmrch+mtZE8C43MymDIii/n5QzgxfygTcjKQY6z3SX8iIswfN4T544bw48/P5OPSQ6wqruSTssNs3FPFrgN1lFTUUlJRG+tQ+5FU99HelrOV1AgW4vpuhTcgIkuBpQB5eT3rMZOcmMDQjLYpKdp/b0nA/RJwX9fnuvEe9fntzu3h+3UScqef3X5/15/XPtaj3yPRI6QleUhL9rTbZqYmkpOZQk5WCsOyUhkxMJUxg9P6f5XEMSzRk8AJY4dwwti2KTpqG5opO1RH2cE6dlfVcehIE4frmqhyHw3NXhrdR0OLl4amFhpbvKg6pWcFvKru67Z9qk5p1+u330RWJJNFGTDG73Wuuy/QOaUikggMxGnoDuVaVPUx4DFw2ix6EuS504Zz7rRze3KpMSaIjJREJg/PsnaofiCSt3mrgUkiki8iycASYEWHc1YA17jPvwS8oc4twgpgiYikiEg+MAmw9RqNMSZGIlaycNsgbgZW4nSdXaaqG0XkHqBAVVcATwC/dxuwD+AkFNzznsFpDG8GbuqqJ5QxxpjIskF5xhhzDAu166y1NhpjjAnKkoUxxpigLFkYY4wJypKFMcaYoCxZGGOMCarf9IYSkXLg0168RTZQEaZwwsni6h6Lq3ssru7pj3GNVdWcYCf1m2TRWyJSEEr3sWizuLrH4uoei6t7juW4rBrKGGNMUJYsjDHGBGXJos1jsQ6gExZX91hc3WNxdc8xG5e1WRhjjAnKShbGGGOCsmThEpHZIvK+iKwTkQIRWRDrmHxE5BYR2SIiG0XkZ7GOpyMR+baIqIhkxzoWABG5z/33Wi8ifxGRQTGMZZGIbBWRIhG5I1Zx+BORMSLypohscn+nbot1TP5ExCMiH4nIS7GOxZ+IDBKR59zfrc0icnKsYwIQkW+6/4+fiMifRaT7C3GHwJJFm58B/1dVZwN3ua9jTkTOwlmT/HhVnQ7cH+OQ2hGRMcBCYGesY/HzOjBDVWcB24A7YxGEiHiAR4DzgWnA5SIyLRaxdNAMfFtVpwEnATfFSVw+twGbYx1EAA8Cf1fVKcDxxEGMIjIauBWYp6ozcJaDWBKJz7Jk0UaBAe7zgcDuGMbi7+vAvaraAKCq+2McT0cPAP+F8+8XF1T1NVVtdl++j7PSYiwsAIpUtURVG4HlOIk/plR1j6qudZ9X43zpRWSN++4SkVzgQuDxWMfiT0QGAqfjrMGDqjaq6qHYRtUqEUhzVxtNJ0LfXZYs2twO3Cciu3Du3mNyNxrAZOA0EflARN4WkfmxDshHRBYDZar6caxj6Q1z7QcAAAM9SURBVML1wKsx+uzRwC6/16XEyZeyj4iMA+YAH8Q2kla/wLn58MY6kA7ygXLgt24V2eMikhHroFS1DOf7aiewB6hS1dci8VmRXIM77ojIP4ARAQ59Fzgb+KaqPi8il+LcQZwTB3ElAkNwqgvmA8+IyHiNUje2ILF9B6cKKuq6iktV/+ae812cKpc/RjO2vkJEMoHngdtV9XAcxHMRsF9V14jImbGOp4NEYC5wi6p+ICIPAncA349lUCIyGKe0mg8cAp4VkStV9Q/h/qxjKlmoaqdf/iLyO5y6UoBniWIxOEhcXwdecJPDhyLixZkHpjyWsYnITJxf0I9FBJyqnrUiskBV98YqLr/4rgUuAs6OVmINoAwY4/c6190XcyKShJMo/qiqL8Q6HtepwMUicgGQCgwQkT+o6pUxjgucUmGpqvpKYM/hJItYOwfYrqrlACLyAnAKEPZkYdVQbXYDZ7jPPwsUxjAWf38FzgIQkclAMnEwkZmqblDVYao6TlXH4fwxzY1GoghGRBbhVGVcrKpHYhjKamCSiOSLSDJOw+OKGMYDgDjZ/Qlgs6r+PNbx+Kjqnaqa6/4+LQHeiJNEgft7vUtEjnN3nQ1simFIPjuBk0Qk3f1/PZsINbwfUyWLIL4KPOg2EtUDS2Mcj88yYJmIfAI0AtfE8E65r3gYSAFed0s976vqjdEOQlWbReRmYCVOL5Vlqrox2nEEcCpwFbBBRNa5+76jqq/EMKa+4Bbgj27iLwGui3E8uFVizwFrcapcPyJCo7ltBLcxxpigrBrKGGNMUJYsjDHGBGXJwhhjTFCWLIwx5v+3d8c2DQUxAIbtDhSoUqWkT0tDySBMwBhI2SQDMAAlDSMhZIp7vcWTTiek75vA3a97UWxaYgFASywAaIkFAC2xgEky83G7qXGTmYft5sB59Vywhz/lwUSZ+RZjz9FtjN1Cl8UjwS5iARNtqyG+YqyQeaqqn8UjwS4+Q8Fcx4i4i4j7GC8M+Je8LGCizHyPcSHvISJOVfW6eCTYxdZZmCQzXyLiu6qu2z3uz8x8rqqP1bPBX3lZANDymwUALbEAoCUWALTEAoCWWADQEgsAWmIBQEssAGj9Ap4aOVWqJGJBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "normal_dist = pm.Normal.dist(mu=0, sigma=1)\n",
    "t_dist = pm.StudentT.dist(mu=0, lam=1, nu=1)\n",
    "x_eval = np.linspace(-8, 8, 300)\n",
    "plt.plot(x_eval, theano.tensor.exp(normal_dist.logp(x_eval)).eval(), label='Normal', lw=2.)\n",
    "plt.plot(x_eval, theano.tensor.exp(t_dist.logp(x_eval)).eval(), label='Student T', lw=2.)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('Probability density')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the probability of values far away from the mean (0 in this case) are much more likely under the `T` distribution than under the Normal distribution.\n",
    "\n",
    "To define the usage of a T distribution in `PyMC3` we can pass a family object -- `T` -- that specifies that our data is Student T-distributed (see `glm.families` for more choices). Note that this is the same syntax as `R` and `statsmodels` use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [lam, x, Intercept]\n",
      "Sampling 2 chains: 100%|██████████| 5000/5000 [00:03<00:00, 1371.34draws/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAE/CAYAAADmL9yLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xlc3NW9//HXYQaYAYZ9S4BAQhJIQiB7YjRm07jlumtd4lKrVa9eb3/em2vb21a93WuvXpdeY1uv9sa1xtbr1d7WVk2MJsZsJpo97AMDzLAPMMDMnN8fwAgEwjaEgXyej0ceAeY73++ZL8P3PWf5nqO01gghhBCBIGisCyCEEEJ0kVASQggRMCSUhBBCBAwJJSGEEAFDQkkIIUTAkFASQggRMCSUhOhGKeVUSk0b63IMhlJqlVLK2u37Q0qpVcPYzwql1DG/Fk6IYZJQEmNCKVWklGrpDIFKpdSLSqmIEewvQymllVLGkZRLax2htS4YyT7GitZ6jtZ660DbdZ6n6d2et11rnTWqhRNikCSUxFj6O611BLAAWAR8b6wKMtIwG+vnCzFRSCiJMae1LgP+D8gBUEpNVkq9rZSqUUqdVErd1bWtUmqJUmqPUqqhs4b1eOdDH3X+X9dZ+zqnc/s7lFJHlFK1Sqm/KKXSu+1LK6XuU0qdAE50+9n0zq+jlFL/rZSyK6WKlVLfU0oFdT52u1LqE6XUE0qpauCR3q9LKfWIUmqLUup1pVSjUmqfUiqv2+NFSqmHlFIHgSallLHztb/ZecxCpdQD3bY3d9Yoa5VSh4HFvY5XpJS6oPNrg1Lqu0qp/M5j71VKpSmlus7Tgc7z9LXuzYCd5dnSa79PKqWe6nZOnldK2ZRSZUqpHymlDAP/loUYHAklMeaUUmnApcD+zh+9BliBycC1wE+UUms6H3sSeFJrHQlkAr/v/Pn5nf9HdzbB7VRKXQF8F7gaSAC2A6/2OvyVwFJgdh9FexqIAqYBK4Fbga93e3wpUAAkAT/u5+VdAbwBxAKvAG8ppYK7PX4jcBkQDXiB/wUOACnAWuBbSqmLOrd9uPM1ZwIXAbf1c0yABzv3fSkQCdwBNGutu85TXud5er3X814DLlVKWaAj3IDrO8sO8CLgBqYD84F1wJ2nKYcQQ6O1ln/y74z/A4oAJ1AHFAP/CZiBNMADWLpt+1Pgxc6vPwIeBeJ77S8D0ICx28/+D/hGt++DgGYgvfN7DazptR9NxwXXALQBs7s9djewtfPr24GSAV7jI8CnvY5vA1Z0Owd3dHt8ae99At8BXuj8ugC4uNtj3wSsvc7pBZ1fHwOu6KdcGpje7ftVvfbzMXBr59cXAvmdXycBrYC527Y3Ah+O9ftJ/k2cf1JTEmPpSq11tNY6XWv991rrFjpqRzVa68Zu2xXTUXMA+AYwEziqlNqtlFp/mv2nA08qpeqUUnVADaC67QugtJ/nxgPBncfuqxyne253vm201l6+qgH2tY90YHJXeTvL/F06woDO53XfvnvZeksD8gdRvr68QkfYANzEV7WkdDrOia1b+Z4DEod5HCFOIZ2rItCUA7FKKUu3YJoClAForU8AN3b27VwNbFFKxdHx6b+3UuDHWuuXT3O8/qbJdwDtdFyID/cuxwDP7S6t64vOMqfS8Rr72kcpUKi1ntHPvmyd+zvUrTz9KaWjme/LQZSxtzeAf1dKpQJXAed022crHbVU9zD2K8SApKYkAorWuhTYAfxUKWVSSuXSUTt6CUAptUEpldBZ66jrfJoXsHf+3/0eo03Ad5RSczqfG6WUum6Q5fDQ0V/1Y6WUpXOAxINd5RiChUqpqztH132Ljov6p/1s+xnQ2DnYwNw5WCFHKdU1oOH3na8npjMw/uE0x/0t8EOl1AzVIbczvAEq6XmeetBa24GtwAt0hOSRzp/bgPfoCKxIpVSQUipTKbVyMCdCiMGQUBKB6EY6+ojKgT8CD2ut/9b52MXAIaWUk45BDzdorVu01s10DDb4pLNpaZnW+o/Az4HXlFINdNQaLhlCOf4BaKKjL+djOpqx/muIr+V/gK8BtcAtwNVa6/a+NuwMwvXAPKCQjtrab+kYbAEdfWnFnY+9B2w+zXEfpyPE3gMagOfp6LODjr6u33Wep+v7ef4rwAV81XTX5VYghI7aYy2wBZh0mnIIMSRKa1nkT4jRoJR6hI4BBRvGuixCjBdSUxJCCBEwJJSEEEIEDGm+E0IIETCkpiSEECJgSCgJIYQIGKNy82x8fLzOyMgYjV0LIYQYh/bu3evQWicMtN2ohFJGRgZ79uwZjV0LIYQYh5RSp5sWy0ea74QQQgQMCSUhhBABQ0JJCCFEwJBZwsVZqb29HavVisvlGuuiCDGhmEwmUlNTCQ4OHnjjPkgoibOS1WrFYrGQkZGBUmqsiyPEhKC1prq6GqvVytSpU4e1D2m+E2cll8tFXFycBJIQfqSUIi4ubkQtEBJK4qwlgSSE/43070pCaZzYtC2fHfmOHj/bke9g07bhrngtJoK33nqLw4cPD7xhL2+//TY/+9nPRqFEw7N161bWr+9Y2X6gstXV1fGf//mfvu/Ly8u59tprR72MY+XOO+8c1u+4txdffJH7778fgE2bNvHf//3fI97naJBQGidyU6O4/5X9vmDake/g/lf2k5saNcAzxUQ2nFByu91cfvnlfPvb3x7Sc4ZKa43X6x3y8wYqW+9Qmjx5Mlu2bBnycU5nOK93NPYB8Nvf/pbZs2f7ZV9d7rnnHm699Va/7tNfJJTGieWZ8Txz03zuf2U/j793jPtf2c8zN81neWb8WBdNDENRURHZ2dncfPPNzJo1i2uvvZbm5mYA3n//febPn8/cuXO54447aG1tBeDb3/42s2fPJjc3l3/+539mx44dvP3222zcuJF58+aRn59Pfn4+F198MQsXLmTFihUcPXoUgNtvv5177rmHpUuX8i//8i89PjUXFRWxZs0acnNzWbt2LSUlJX0+p7sXX3yRK664glWrVjFjxgweffRR376ysrK49dZbycnJobS0lPfee49zzjmHBQsWcN111+F0OgH485//THZ2NgsWLOAPf/hDj313la2yspKrrrqKvLw88vLy2LFjB9/+9rfJz89n3rx5bNy4kaKiInJycgBYtmwZhw4d8u1r1apV7Nmzh6amJu644w6WLFnC/Pnz+Z//+Z9Tfidbt25lxYoVXH755b4QeOmll1iyZAnz5s3j7rvvxuPxAPD8888zc+ZMlixZwl133eUrb+9z1t9xDx065Ntvbm4uJ06coKmpicsuu4y8vDxycnJ4/fXXe7wGgFdffZW5c+eSk5PDQw895Ct7REQE//qv/0peXh7Lli2jsrLytO+/Rx55hF/+8pe+/T/00EMsWbKEmTNnsn37dgA8Hg8bN25k8eLF5Obm8txzzwFgs9k4//zzmTdvHjk5Ob7t/UZr7fd/Cxcu1GJ0/Ptfjur0h97R//6Xo2NdlHHt8OHDY3r8wsJCDeiPP/5Ya63117/+df3YY4/plpYWnZqaqo8dO6a11vqWW27RTzzxhHY4HHrmzJna6/VqrbWura3VWmt922236TfeeMO33zVr1ujjx49rrbX+9NNP9erVq33bXXbZZdrtdmuttX7hhRf0fffdp7XWev369frFF1/UWmv9/PPP6yuuuKLP53T3wgsv6OTkZO1wOHRzc7OeM2eO3r17ty4sLNRKKb1z506ttdZ2u12vWLFCO51OrbXWP/vZz/Sjjz7qe53Hjx/XXq9XX3fddfqyyy47pWzXX3+9fuKJJ7TWWrvdbl1XV6cLCwv1nDlzepzLru8ff/xx/YMf/EBrrXV5ebmeOXOm1lrr73znO3rz5s2+czdjxgxfmbp8+OGHOiwsTBcUFGitO94j69ev121tbVprre+99179u9/9TpeVlen09HRdXV2t29ra9Hnnnecrb+9z1t9x77//fv3SSy9prbVubW3Vzc3NesuWLfrOO+/0laeurk5rrfXKlSv17t27dVlZmU5LS9NVVVW6vb1dr169Wv/xj3/UWmsN6LfffltrrfXGjRv1D3/4wz5/Z13lfPjhh/Vjjz3m2/+DDz6otdb63Xff1WvXrtVaa/3cc8/59uNyufTChQt1QUGB/uUvf6l/9KMf+X4nDQ0Npxyrr78vYI8eRH7IkPBxZEe+g5d2lfDAmum8tKuEZZlxUlPyA4/HM+Any6FKSkrCYDCcdpu0tDTOPfdcADZs2MBTTz3FhRdeyNSpU5k5cyYAt912G7/61a+4//77MZlMfOMb32D9+vW+/pfunE4nO3bs4LrrrvP9rKuWBXDdddf1WaadO3f6aiq33HJLj1pRf88BuPDCC4mLiwPg6quv5uOPP+bKK68kPT2dZcuWAfDpp59y+PBh3+tsa2vjnHPO4ejRo0ydOpUZM2b4Xv+vf/3rU47xwQcf+Po+DAYDUVFR1NbW9lkegOuvv55169bx6KOP8vvf/97X1/Tee+/x9ttv+2oHLpeLkpISZs2a1eP5S5Ys8Q1lfv/999m7dy+LFy8GoKWlhcTERD777DNWrlxJbGys7xwdP368z3PW33HPOeccfvzjH2O1Wrn66quZMWMGc+fO5Z/+6Z946KGHWL9+PStWrOhRtt27d7Nq1SoSEjrmNL355pv56KOPuPLKKwkJCfG9JxYuXMhf//rXfs9RX66++mrfc4uKinxlP3jwoK9ptL6+nhMnTrB48WLuuOMO2tvbufLKK5k3b96QjjUQCaVxoqsPqavJbllmnDThjXO9RymdbtSS0Wjks88+4/3332fLli0888wzfPDBBz228Xq9REdH8/nnn/e5j/Dw8CGX8XTP6a/83Z+jtebCCy/k1Vdf7bFtf2UcqZSUFOLi4jh48CCvv/46mzZt8pXjzTffJCsr67TP71322267jZ/+9Kc9tnnrrbeGtI++jjtr1iyWLl3Ku+++y6WXXspzzz3HmjVr2LdvH3/605/43ve+x9q1a/nBD34wqNcdHBzsO/8Gg2HI/VmhoaGnPFdrzdNPP81FF110yvYfffQR7777LrfffjsPPvigX/unpE9pnDhore8RQF19TAet9WNcsvHPYDAwefJkv/4bqJYEUFJSws6dOwF45ZVXOO+888jKyqKoqIiTJ08CsHnzZlauXInT6aS+vp5LL72UJ554ggMHDgBgsVhobGwEIDIykqlTp/LGG28AHReVru1OZ/ny5bz22msAvPzyy6d8Qu/PX//6V2pqamhpaeGtt97y1Ya6W7ZsGZ988onv9TQ1NXH8+HGys7MpKioiP79j9Gjv0Oqydu1ann32WaCjRltfX9/jNffla1/7Gr/4xS+or68nNzcXgIsuuoinn34a3bnS9v79+wd8fWvXrmXLli1UVVUBUFNTQ3FxMYsXL2bbtm3U1tbidrt58803+91Hf8ctKChg2rRpPPDAA1xxxRUcPHiQ8vJywsLC2LBhAxs3bmTfvn099rVkyRK2bduGw+HA4/Hw6quvsnLlygFfx3BddNFFPPvss7S3twNw/PhxmpqaKC4uJikpibvuuos777zzlHKOlITSOHHPysxTakTLM+O5Z2XmGJVIjFRWVha/+tWvmDVrFrW1tdx7772YTCZeeOEFrrvuOubOnUtQUBD33HMPjY2NrF+/ntzcXM477zwef/xxAG644QYee+wx5s+fT35+Pi+//DLPP/88eXl5zJkzp88O/d6efvppXnjhBXJzc9m8eTNPPvnkoMq/ZMkSrrnmGnJzc7nmmmtYtGjRKdskJCTw4osvcuONN5Kbm+trujOZTPz617/msssuY8GCBSQmJvZ5jCeffJIPP/yQuXPnsnDhQg4fPkxcXBznnnsuOTk5bNy48ZTnXHvttbz22mtcf/31vp99//vfp729ndzcXObMmcP3v//9AV/f7Nmz+dGPfsS6devIzc3lwgsvxGazkZKSwne/+12WLFnCueeeS0ZGBlFRfY+C7e+4v//978nJyWHevHl8+eWX3HrrrXzxxRe+wQ+PPvoo3/ve93rsa9KkSfzsZz9j9erV5OXlsXDhQq644ooBX8dw3XnnncyePZsFCxaQk5PD3XffjdvtZuvWreTl5TF//nxef/11/vEf/9Gvx1VdCe5PixYt0rKekghkR44cOaU/4UwqKipi/fr1fPnll2NWhpF48cUX2bNnD88888xYF2VMOJ1OIiIicLvdXHXVVdxxxx1cddVVY12sgNHX35dSaq/W+tRPLr1ITUkIIYbokUce8Q2Jnjp1KldeeeVYF2nCkJqSOCuNdU1JiIlMakpCCCEmBAklIYQQAUNCSQghRMCQUBJCCBEwJJSEGAO9Z7o+Wyxfvtwv+7n99tt909/4a2kHERgklIQYA6cLJX8teeDP/XXNjj1SO3bs8Mt+uhuNpR3E2JFQEmIM9F5+ofeyCd2XYwD45S9/ySOPPALQ7/IU3T3yyCPccsstnHvuudxyyy39LkPg9Xr5+7//e7Kzs7nwwgu59NJLfTWQjIwMHnroIRYsWMAbb7zR73HfeOMNcnJyyMvL4/zzzwf6XpoBOpZYgI4pkDZu3EhOTg5z5871LdOwdetWVq1axbXXXutb2mOg21a6L+3Q3xIOdruda665hsWLF7N48WI++eSTof/SxJkxmKnEh/pPlq4QgS4Qlq7ovvxC72UTej/+2GOP6Ycfflhr3f/yFN09/PDDesGCBbq5uVlr3f8yBG+88Ya+5JJLtMfj0TabTUdHR/uWwkhPT9c///nPffvs77g5OTnaarVqrb9aUqOvpRm01jo8PFxrrfWWLVv0BRdcoN1ut66oqNBpaWm6vLxcf/jhhzoyMlKXlpZqj8ejly1bprdv337K6+u+ZEfX0g5a97+Ew4033ujbT3Fxsc7Ozu7vVyP8QJauEGIEtm7tf3bukVq1avA3p3dfNqE/Ay1P0d3ll1+O2WwG+l+G4OOPP+a6664jKCiI5ORkVq9e3WMfX/va1wY87rnnnsvtt9/O9ddf71sCoa+lGbr7+OOPufHGGzEYDCQlJbFy5Up2795NZGQkS5YsITU1FYB58+ZRVFTEeeedd/qT16m/JRz+9re/9eh3amho8E0VJAKLhJIQAaL7kgdGo7HHUuIulwsYeHmK/van+1mG4E9/+tOg9nG6427atIldu3bx7rvvsnDhQvbu3ctNN93U59IMg9G1jAIMfRmG/pZw8Hq9fPrpp5hMpkHvS4wN6VMSYgwMtPxCUlISVVVVVFdX09rayjvvvAMMf3mK/pYhOPfcc3nzzTfxer1UVlaydevWPp9/uuPm5+ezdOlS/u3f/o2EhARKS0v7XJqhuxUrVvD666/j8Xiw2+189NFHLFmyZMDXMVzr1q3j6aef9n0/Wus5iZGTmpI46w2lic1fui+/cMkll3DZZZf1eDw4OJgf/OAHLFmyhJSUFLKzs32Pvfzyy9x777386Ec/or29nRtuuIG8vLzTHu/OO++kqKiIBQsWoLUmISGBt956i2uuuYb333+f2bNnk5aWxoIFC/pdhqG/427cuJETJ06gtWbt2rXk5eXx85//nM2bNxMcHExycjLf/e53e+zrqquuYufOneTl5aGU4he/+AXJycl9Dtrwh6eeeor77ruP3Nxc3G43559/vm8BQBFYZEJWcVaSCVm/0tW3Ul1dzZIlS/jkk09ITk4e62KJcWwkE7JKTUmIs9z69eupq6ujra2N73//+xJIYkxJKAlxluuvH0mIsSADHYQQQgQMCSVx1hqN/lQhznYj/bsaVCgppf6fUuqQUupLpdSrSikZ7C/GNZPJRHV1tQSTEH6ktaa6unpE94MN2KeklEoBHgBma61blFK/B24AXhz2UceJTdvyyU2NYnlmvO9nO/IdHLTWc8/KzDEsmRip1NRUrFYrdrt9rIsixIRiMpl8M3IMx2AHOhgBs1KqHQgDyod9xHEkNzWK+1/ZzzM3zWd5Zjw78h2+78X4FhwcPOCUPkKIM2/AUNJalymlfgmUAC3Ae1rr90a9ZAFgeWY8z9w0n/tf2c+GpVN4aVeJL6CEEEL434B9SkqpGOAKYCowGQhXSm3oY7tvKqX2KKX2TKQmkeWZ8WxYOoWnPjjJhqVTJJCEEGIUDWagwwVAodbarrVuB/4AnLJ8pNb611rrRVrrRQkJCf4u55jZke/gpV0lPLBmOi/tKmFHvmOsiySEEBPWYEKpBFimlApTHdPvrgWOjG6xAkP3PqQH12X5mvIkmIQQYnQMGEpa613AFmAf8EXnc349yuUKCAet9T36kLr6mA5a68e4ZEIIMTHJhKxCCCFG3WAnZJUZHYQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAGFUpKqWil1Bal1FGl1BGl1DmjXTAhhBBnn8HWlJ4E/qy1zgbygCOjV6SJZ9O2fHbkO3r8bEe+g03b8seoREIIEZgGDCWlVBRwPvA8gNa6TWtdN9oFm0hyU6O4/5X9vmDake/g/lf2k5saNcYlE0KIwGIcxDZTATvwglIqD9gL/KPWumlUSzaBLM+M55mb5nP/K/vZsHQKL+0q4Zmb5rM8M36siyaEEAFlMM13RmAB8KzWej7QBHy790ZKqW8qpfYopfbY7XY/F3P8W54Zz4alU3jqg5NsWDpFAkkIIfowmFCyAlat9a7O77fQEVI9aK1/rbVepLVelJCQ4M8yTgg78h28tKuEB9ZM56VdJaf0MQkhhBhEKGmtK4BSpVRW54/WAodHtVQTTFcf0jM3zefBdVm+pjwJJiGE6Gmwo+/+AXhZKXUQmAf8ZPSKNPEctNb36EPq6mM6aK0f45IJIURgUVprv+900aJFes+ePX7frxBCiPFJKbVXa71ooO1kRgchhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgQMCSUhhBABQ0JJCCFEwJBQEkIIETAklIQQQgSMQYeSUsqglNqvlHpnNAskhBDi7DWUmtI/AkdGqyBCCCHEoEJJKZUKXAb8dnSLI4QQ4mw22JrSfwD/AnhHsSxCCCHOcgOGklJqPVCltd47wHbfVErtUUrtsdvtfiugEEKIs8dgakrnApcrpYqA14A1SqmXem+ktf611nqR1npRQkKCn4sphBDibDBgKGmtv6O1TtVaZwA3AB9orTeMesmEEEKcdeQ+JSGEEAHDOJSNtdZbga2jUhIhhBBnPakpCSGECBgSSkIIEeA2bctnR76jx8925DvYtC1/jEo0eiSUhBAiwOWmRnH/K/t9wbQj38H9r+wnNzVqjEvmf0PqUxJCCHHmLc+M55mb5nP/K/vZsHQKL+0q4Zmb5rM8M36si+Z3UlMSQohxYHlmPBuWTuGpD06yYemUCRlIIKEkhBDjwo58By/tKuGBNdN5aVfJKX1M/jSWfVgSSkKIs9J4GjzQ1Yf0zE3zeXBdlq8pb7SCaSz7sCSUhBBnpfE0eOCgtb5HH1JXH9NBa/2oHK97H9bj7x3zBeKZaDJUWmu/73TRokV6z549ft+vEEL4U1cQTfTBA8P1+HvHeOqDkzywZjoPrssa0b6UUnu11osG2k5qSkKIs9bZMnhgOM5kH1Z3EkpCiLPWWF14A92Z7sPqTkJJCHFWGssLb6A7031Y3UmfkhDirLRpWz65qVE9mux25Ds4aK3nnpWZ4+YY44X0KQkhxGncszLzlD6k5Znxfg2L8TTCL1DINENCCDFKzqbpgfxFakpCCDECA92EKyP8hkZCSQghRmCgJjp/j/Dz10wUgTqjhYSSEEKMwOlmPxiNEX7+6qcK1P4uGX0nhBB+0NfsB6M1+s5fM1GcyRktZPSdEGJCCsRmp/6a6EZrhF9//VRDPTeB2N8loSTEBBWIF29/CLRmp7G4Cbe/EBzquelrP1prKisr+eyzz3jrrbfYvXv3qL2OvkjznRATVPeLZe/+jbH+RDzSZq1Amkj1TN8gO9DvdbDnZke+g3ue/4hvLY0ixeTmowMnePH9g1ycHU1OxiTi4mKZPNnMlCma+PgYYmPXjajcg22+k1ASYgIbzsX7TFxk/RGY/pzBejwZzO+n97lpbGzk5MmTNDQ0UFpaSlVVFR98UUx6YjQLs9JISPByyH6SSHMZZmMp0cYCzOYqTKZ2AMLDc1m8+MCIyi2hJIQAhn7xPlM1rJHUdgKpphRIWlpaeOODz/jxWwdYmuTlg73HWDfDQlZKHAkJ0URHO4mKqgeKcbvzaW4+isVSi9HoPe1+lQrl/PObUMow7LINNpRkRgchJrDefQbLMuMGvHifqVkIuneyP7Bm+pADqatMyzLjAqZZ8kxpa2vj5MmTNDU1YbPZKCwspKamhgpnOx8WN/DNCyOYm9HGxQtPgAaMAAAgAElEQVSslFcfYlq8g+goJ0FBX1VCDAYIDR3gOO0hGA1TmTRpGR5PE0Zj5Ci/MgklISaskVy8hxsYQy3fUAMTTj+D9UQLJY/HQ0FBAfX19dTV1XHkyBFqamowGAxMmmQmPLwas7mSjIxipkwpQhnKuTGyxff8yaEwOXbg4zS1hNPmisfhCOdwUQi7jngJ9SSQmZDOmjVrWbPmplF8lT1JKAkxQY3k4j3cwBiskQRmX/1ayzPjx3Ugeb1eSkpKcDgctLS0cOjQISorKwkKgkmTjFgsNXg8BSQmlpGaasNstmMytQ75OE0tMRwrCcHbEovLGc2xQsW2I8HMnzKFBbMyaTKE8/pnR1izaDa7q0PYcOkszl+UPQqvuH/SpySE6OFM9CmNxZIOgbCMhNYam81GWVkZXq+XI0eOUFxcjFIeJk3yEhFRg9N5CLPZTmRkLWazg+Bg95CO4fUG0dqagMMRRm1tJBUVIVRXW/B4JmMyRaHD4/jD3hKWzZzMl65Yvn15Hnmp0ew8WsrPXv4Ll8yMJD3eQkFVA//76VE2XLSM//jRv474tUufkhBiWM5E81jvEOgKjO4/93dgdN3D01fYjgaHw0FRURFaawoLCzl27BjgIjHRRWhoBcUVn5OS7GThwibCwmp69PcMhscTgtMZS01NBA5HOOXlRlyuZLzeZCyWaCwWCwaDgQULsklPTwfA6XRy5MgRTpQ5+PBAIUszGzi+s4q/lZRwoqqBtVkzWZozg6qqKkxuJ5efM4uQFKkpCSHOMr1rY9/5w0HeOWjjuVsW+oIwkKbn6a6hoYH8/Hw8Hg8VFRXs27cPpZxERdViNJbjducTG+skLs6J2Tz0lVvb28Oor4+ipsaC1WqkttZCe/tk3O4YIiIsxMTEEB8fT15eHiaTyTf8u6KiAq+3Y1Sd2+3GZrPR1NREUEQ8W0vbyAxzcbjEwZo5qaxckovFYvEdc+bMmSxcuJDg4OARnZvupKYkxFkmEJqnhqv3iL93Dtp6PO6vWs1IBnA0NTVx8uRJ2tvbqa2tZe/ePXg8lRgMVszmKgwGK3FxTSxd6iQ0tHnIZXPUmXE2RlNfE055eTBOZxwuVyKtrWZSUlKIjo5m8eIcJk2aRF1dHUVFRVRXV9PS0oLD4eC9996jsbGR+vp6IiIiCA8Pp7q6GgCj0ciUKVNwmWJ4ZZeVB65awtevuoiTjUHc/8p+vnVx3+Hc3t6Oy+XqEVijTUJJiAlitJqnzlTY9Q6MrsEP/qzVDGYAR2trKydOnMDlctHU1MS+fbuprz+CUqVYLDWEhdmJjXWycKGT4OC2IR3f61U4nZHU1looLw+hujqcw8Umdh4NYmZqKhcvmUNKSgrr18/C5XJRXFxMU1MTDQ0NtLW1sXfvXpqbm2lpacHtduNyuXy1maCgIKKiosjLyyMuLo65c+cyd+5cgoK+mk1u07Z8Xrnmq99lQoLmhxel8s4HnxBSm3JKeY1GI9OmTRvqaR4RCSUhJojRur/oTPXF9BUY/hyW3ruJcFF6FHc+9Tb/tHoKWfEm9u3bSXn5ZwQHlxMT04jFUkNsbCNzc50YDae/ubQ3t9tIfb0FhyMcmy2E2loLtbWRtLUlkpycyqRJk8jJmc1xRwuvH/+ENati2Hm0FGu9i6AgG3/c/jkWgweLwU1zczMWi4XKRhfVzjYWZSaTkJBATk4Oa9euJSrq9HP+NTU1UVRUhMvlYnEE6Npa9uwp8j0+NzmRixeux2AY/o2x/iR9SkJMMKMx/c5oz6DQ14i/uzfvBSA3JYqDZfXD7l/yer0UFRXxm78eJC3SgMu2H5vtU8zmKsyRDURaaklNbiEysgmlhlbu1tYQamoisNvNVFSEcrQomDZXAvGR00hNTSMzMxNbq4EvjhWxelayr9ajtabQ3sg7nx1j4SQTCZZQalu9fFbcwLmZcRiNwXxS1s5Dt/0d9910BXtKG/odAel2uyktLcXhcKB6vYA391rJnZrElSvmYTabh3zu/En6lIQ4C43W/UWjfTNt7xF/XdbnTuLv8iZz9+a93L15L8/dshCgz5qa1pqysjJsNht2exUHDnxITc1+wsIcxMU1MTumnthIJ+GTXEMuX0OjiYLSYKrtZo4WKpQ7kVCVRnh4KtOnzyA+Ph6DoQlPeCsvf1rEJdFhGI1G3v/0c9759BBL0qPZ1VCM1+v1NbcVVji5+cJlXH/pKmbMmEFiYiI7C6p9gdEV1N4PC3jxg4M8uDyOkNqiHrUc6GhiS0tLIyMj45RQaovJ4P5X9pOe2cTyTPOojzj0B6kpiQljPHf0+8No3l90puea6/277Ko5zU2J4oitgR9elEqsp46CgpMUFu6kqekw4eHVxMU5iYlpIDbWSWho+5CO6fFCba2ZakcYRaUGDhcYMHqTsZgyiY9P51h9ENtPOFg3byqX5KVSVVXlG2TgcDgIDg7GZDJR1dDKJ/kOMhMiKKhp5e5r1vF3KxcRFhZGZmYmYWFhpxy7exNbF601L39axGu7rXzjglwevWnlsJrYAmWeQKkpibPOmb4PJdAMdH/RcEN7LOaa6ypPTU0Nn376KWVlhWSXv0d1/kkuzHBR9nkrrjgnyclOUlM9Q9q32x1EdbUZmy2EqioTtbWRuFxJ1LROZmeJJm9KLEedwdy+YiZxwe2UlZXx+bEC9hTVMCM5hve376StMpGU2HDa29uZNGkSK1euJDU1FaPRSFxcHH845uKZbYV8p7MJtauJ7fDhw6eURylFWFgY06ZN8zWxQcd531Zby8ab5/LSrhIuLaod1vk+E1NG+ZPUlMSEEiifCgPRcGtSpwszwC+107KyMrZv347Nlk9t7QFCQmyYzVVERdWTkNCExdI05JtLXS4jlZUhVFSE4nCE0dAQg8eThlLJmEzhZGZmEhsbi9Pp5IsvvqCqqoqi+naO2NvITYkkNzUak8mEvcnN9grFN9fNY/WCbI47XDz23nEev+MCLl3c0WentcbhcFBaWsr+khoe+/MxLslJ4v++rGTjxVkszIgnLS2N2NjYU5rY+uLPWm+g/E3I0hViQhjOp/uzdZ2dwfD3BWooF8/29naOHDnC3r17qKg4Snv7SUJDKzGZKomOric+vpnw8KHf39PYGIzNFkxlZUetp6Ehhvb2yXi90cycmUVKSgo1NTUcOHCApqYmLBYLERERuN1utNaEhoaSnZ1NudPLy7vLuHhxFlsLm9l4cRZLpyfzQbliUWYSeclmiouLaWlp4YC1jhMVjVyzMNVXjvj4eEpaTXzr91+MOEz81RQdSAs9SiiJCWGof1SB8qnQn/zdV+bv0O5+zjd/Wsy/rUsh1lPHvn17cDoLgGJCQipwB9mYFN9IYkITJtNQ7++B2toQysuDsdvN2OwRlFZGEBWcxuTJWUyZMoXq6mrffT1RUVHExsbS3NyMx+MhJiaG5cuXo7WmoaGB2NhYEhMTCQoKQmtNaYuRf3vfxvfXJJMe7uWAtY7H/nyMjRdnkZcaDUB4eDjp6el99gl1CbR+zUAqj4SSmDCGsrxzoHwq9KdAbMppamqioKAAm83GF1/sY2/RF1Q7j7NwWgsz01oJC7MTE+MkOHio/T0Kuz0Emy0Yuz0MhyOc5uYEwsKyMJujaW1tpbW1FafTSWRkJImJiTidTrxeL1OmTCEnJ4fIyEisVishISGkpKQQHByM1hqDwUBERARNTU14PD3L9eZeK1mTIlm/PJe4uDiUUmfVIJkzQUJJTCiD+XQfSJ8K/c0fYTLUcPN4PBQXF1NZWYnVaqWs7CRudyFBQaWYTFVER9djMlUSFdWEwTDU/h5FsdWIvTIUe5WJurpIjMZMampCCQ0N94VHWFgYycnJtLS0EBoaSnp6OgkJCUyaNImGhgZqa2uZPHky4eHhaK1pbW2ltraW+Ph4QkJCgI7+ni7x8fG+AQnizJJQEhNG9wvyb7YX8uC6Gdy1YvRmkw5UI2126y+0d3xZyOrUIOx2O5WVlTQ12WhtPYbBUIbFUkN4uAOTqYrIbovHDVZjowGbLZiyMiOVlaG4XEkcKjLyZamRnMx05kyOJCQkhOTkZIKCgoiOjiYmJgaj0cjkyZMxmUzYbDbCw8OJjY3F7XZjt9upq6sjNjaWpKQk38CBrlFsAzWxibEhQ8LFhND707zFbOQn7x4F4K4VmT0eH696h8WmbfkYgjrum+kK2t9sz+c32wuHfVOsy+ViRbyLWpuN13Zb+fPBUiZH1pBkthJpruJIox1lrCA2zkF6+tDu7wGorjZSXm6ktDSI8koTVTWxOGoiiYpMITMlsTN4EjnZaKS0pparrklkR34Nc87N4/zcTFwuF06nk7i4OBobG6mqqqK0tBSDwUBKSgoREREAmEwmli9f7mtiExOPhJIIaL3vvemqIT3+3gkaW9wTYjBD7/urDEHwk3eP8t3LOtax+c32fN/3d63I7Pc+Ia/Xi9VqpaKigsrKSmpqalDKS3t7CR5PPuHh1ShVgslUxY0rqwkNHVp/j8cDFRUGrFYDZWVGGhqiaWlJpKEhisTEdOLj40lOTmblytm4XC4+L6xg885i1uRmsXpBFn86bOdvHx7lpjXzOCfFxPQkC7/68+c01NeSl5FESkoKRqOR6dOns3r1amliO0tJ850Yl0arKWusmgF79xndu2oaz24t6LPJsq6ujre27+dAfjmLExUejwetXUApBoMVg6EMrYsxmSqxWBowGoc2mWhrK1itQb4mt8bGWJqb4/F6JzN7di7Z2dlkZmZSV1eH3W4HwGAwkJiYSExMDG63m4qKChoaGjhR1cjmncUsz4zjw2MOrjtvFndfupSMjAzCw8PPmqZXIc13YgLzx/xugTb7Q++77m9bmkb+8eP8YvM7XDw1hNDCKn5X8DFBQc2Eh9uJUlaWJh8HiomIcGCxOIc8mWhDA5SUBlFUYsDuMONtS6KlJZHc3AtYtWoNeXlhOJ1OysrKfKPVoqOjSUpKorGxkebmZoKCgvB6vXg8Hmw2Gy9/eIBpSVHccvkFvrnYnK9/zh/3l7HxvlM/QCzPjB9RLfdsXVZ9IpNQEuOKv6a8Ga1lHoZCa01FRQVlZWV8fKiIp97ezbKpMTz13FYKPwunsKmcB86poqnlOMpYS3x0HWFhQ59M1OFQWK0GKio6ZrRubU0mK2sdSbPW8fQXx7hwZhwfFR1mw5JJzEiIwGQy0dbWRmVlJUFBQVgsFurq6tBaU1dXR11dHVFRUcyaNYsZM2b0WJ10eufv5wKvhalK8Zvt+by1v4yr5qf4dYLYLmPx4WK4x5QwG5wBm++UUmnAfwNJgAZ+rbV+8nTPkeY7MVoC/UbSvjidTgoKCnA6nVitVlpaWvjwmJ0pcSaWzQqlps3Kx4e2kZNURVRYJfGxDZhN7iEdw+OBysogysoMVFaacDrjyMg4n3nzLicmpmPWgcrKSmprawE4UdnI5l2l3LpsCnPS4qhsD+E/3tnHhsWTmJFowWQyMWnSJMxmM6mpqT1GuQ2k6yK9cmY8b+0v9/WFjdZ9Y2MxOnM4Q/Qn6n10g+W3IeFKqUnAJK31PqWUBdgLXKm1PnVmwU4SSmI8GM6Fpb9Q3F9UzcUZRmpqarDb7djtdpRSKKUwm41MmuShvf0EdXUHcLmOYwi2ERVZT0jw0Pp029qgrCwImy2YhoYYzOZspkxZQUrKEmJikqitraWyshKv14vT6aS1tZWYmBiCg4NJTk4GoKKigg+OVpI7I53z5mQAEBISgl1FcaLG7ZeLd1fYXzV/Mk987asaxGjVDLofb9txxxm58A/nA81EnHFksPzWp6S1tgG2zq8blVJHgBSg31ASgU2aEYbXDKi1Js3cxh2/fJ0Hzk8jhiYOl9eyeUcxty3P4NkvjcxNa2SypYyIiC/xegsxhlQQZWmkvXOUdXT04MrndEJ5uRG7PQyl0gkJmUFc3AIyMpYyfXooZnMZjY2N1NTUYLN5aG8vx2arJjY2lhkzZlBVVUVUVBTJyck9RrHFxMSwevVqNvRaAqHjPRHOPXNH/p7o3ee3I9/RY+Zyf1+Eex/v3lXTRv3CP9x+zfE2Y/dYGFKfklIqA5gP7BqNwogzI9A6+YdrJOF6umUe5k8OJz8/H6fTSWVlJXV1db7nxcbG8C8XxvLm7nc5L6MGT9BxfvK1OqIjHISFfTWZaFzc4F5DTY2iyGqkri6aCFMmSmUQH7+QqVMXERNjx+kspbGxEY8nCKczkpKSMt8cbCEhIXg8HoKCgnrsUynFOeecQ0JCwuAKgf/eE2d6mYv+jrdyZsKoXfhH8hpHaxHGiWTQQ8KVUhHANuDHWus/9PH4N4FvAkyZMmVhcXGxP8sp/GwiNCOMpI3e6/VSWlqK3W6nsbGxxwizjv6UZJSqxOk8TGvrcaAEKMZsriI0dOiTidrtBqqrI3A6Y/F4UoEpuMPm89vtDmaFNfG5tY7LcpLJSo0nISGBkJAQ0tLSyM7ORilFaWkpXu9XQ7u11pjNZjIzM3uswTMS/nhPnOlaeF/H+832fB5/7wR3rZg6Ku9tf61LJX1K/Ww3mFBSSgUD7wB/0Vo/PtD20qd05gX6Eg+jdbEa6EJaW1vrW9GzrKyMpqYmAIKCgkhOTsZiMdPQcISmpkMEBZUSFGQFSjCb7RiNQ7u5tL0dymxGKu3hnCwLR5FKetxsQkKm43J5faFiNBqp1yZ+v7+K7964mtsvX8OuolruefYvPLhysm9WaugIn4SEBKZMmXJKjWg0jPdlPwL5wj9afwPjpTnenwMdFPA7oEZr/a3BHFxC6czz1xIPox0eo3Gx+MU7X/DkW59wbW4cl2ZFUllZyQdHq0iLNbNoRsfIsba2evaf2EF983HmpTdhMFhRqhSTqXrIi8c1NysqKkNpqI+kpCKcz4vM4EnhqD2GRZmT2V1UQ25KFIerPdy5bh5XrV7K7NmzMRqNNDc3k5+fzys78pmeFEFuShRKKYKCgqj0RlDcZBizC8lEqD2Plwt0XyZ6DcyfoXQesB34AuhqP/iu1vpP/T1HQmls+GOJB2DU3uAjuehprSkvL8dms+FyubBarbS3t3O8qpGXPrVyydJZvH+8hnuWmZgcXkpNSz6ldSdYkNFEeGgZZnPjkMtbXx+E3W6mojKMLwoMWEKmUlodTULsdKYnxwIQFhbGJ9VmPrMHce3qBXxS5OTpG+cxI1Lzzs4v+MX/HWXjRTPJS4vxbZ+ZmUloaOiQyzOaxsuFbSIbye9gPHygkFnCz1L+WOJhNN/gA5WvoaGBgoICWltbsdvtVFdX+x5LSEjonLCzgerqo9gajnPYdogVWS4soWUYg21EhA99MlGHo2M6Hbs9jLq6SLzeKZTaI/nEGsr8qYl8aW/nX264gPNzM0lPTyc2tiOQ3G43b27by8Nv7OHiOUm89Xk5Ny+bwlXzU0lMTCQ1NZVPC2sm9Kf0QDfeXtdI/vYCvelVphkaJn+9icfij2GwI3v6On73obqjNWy1q3z3nZ/Ob/+0k6jmUjKjjZSXl/uWpg4PDyclJQWv10tbWzNGYznBweW0t+dTXV2I01lJdHQDSUkekpJg3ozBHdvthsrKr8KnqSme1tZkIBVrazRzZqSzasksEhMTSU5OpqglBPPWAnbkV3P/NZNYlWGgra2NgoICCgoKAPjS1sjjO+v4r41fY3lmPNd3XlAui8lgypSvmkK7n+9AvSAO9J4Yr8bb7AvD/dubSKP6JJR66f4mPmitxxAEz24t8L2JB/vGPNPDrv05FNcfb3CtNVVVVb5mtg/3H+e5v33JhqVTmFxbxb2Lo/nF/x3l7qWJJIY04/EU0tJyBKezmIaGSqKiGoiOdhIbO/TF4yoqjNhsoVRUmTlcGoZ2p3C4OplbLjmPxTPTWL48ldjYWKZOnUpwcDA78h3c9/I+zouZDLj4y56j/PRPRwC4Km8ym9/fz/J7L2Fl3uQex9qzLZ9N35jb57Dy5Znxo/IeGMwyF4EafGNhuNNJjdVtE8P52zvTw/BHmzTf9cFf06ScyXZef32yG2q7dlfHvcvlor6+noqKCjweD0op4uLiSExMpKWlhf/6y26iVSVxocUEB5djMlURFlFLXJyT2OiWIU8m2tioKC8PpqIilJoaC62tSURFzcViyWTSpBSaQmN5dmcVT999ERfkZvhex39cl0OSaqChoQGllG9V0oNl9TzxiZ2vXzCfF3d23M7w3C0LR9y/4u/3QO+y9F7WQvqC+jYeZl8Ybp/SeGmilD6lEeo9bclw35iB3s7bW19v8O3HK9m27wgXZYbjdrt9Q6uDgoIIDQ0lJSUFg8HAsWPHOHbsKK2tZZjNdkJDK4iKqusInlgnZvPQJxOtru5YQqGqykRTUzxebxrh4bMJC5vErFmzSU9PJz29Yy2f7nOzPfGnz4lx15AzObJH8BQ4Wnjo+pVERkaecqyu39W5mXHc16vpJJDm1zvdMheB2sk9lsZLP814CZfhkj6lEehdhR7u3eHjrZ23urqaJZZ6PDU1/PnEHqqqqoCOe3qWJiXR2trK0aNHOXToEM3NjYSF1RMcbOP48VpiY50kJrawalUdISFDG2zg9XYsHldebqS6OgKXKxG3O4WoqLnExKQwf/58pk2bRlpaGoZu0+NorSkrK6O4uJiioqIeoXTZ9EimTp3TYwbrxYv7L0Pv31WX7heKrt/dUC4Uo/Ee6N3vcNeKTBpb3DJ1TR/G0+wLE7Vfb6gklHrpb/ntoU69P9g/hjP96cjlclFQUEBzc7NvaHVbW8cMBREREYSEhFBYWMjx48dxuVwo1Y7BUI7FUkt8fBPJyW2sXl1HeHgtBsPQF48rKzdQYQumvj6K5uYEtJ5CSWMac7Jncf3F5zNr1ixMJlOPc9DW1sbJkyfZt29fj+BRSpGSksLChQtHtDT26X5XI+lbGK22/t4XS4vZOK4+/JxJp5tO6nTnaKL104wn0nzXS/eQ6Hpj3rtqGh7vqZ2fg91Pl77CZjTuD/F6vZSVlflmi+6+ZIHb7cbr9VJbW4vNZqO+vp6goCBaWhxYLDUkJ7eRlOQiNrYRs9lOeHjjkG8ubWpSlFqDKCwx0OyMpb01ifq2NLaVxnHP+vO4+4a/Iykp6ZRz8JPLppKkGtlXXM1jfz7GxouzyEuNJiQkhGnTphERETGs8zGQ0RoiPxofOKRP6cyY6E1pY0H6lPzgTA0PH+5Fr76+nsLCQtxut2+F0La2NhwOB263G7fbjcvl8n3f2NiA2ewiNdVDYmIL4eHVREbWYbFUYzY3D3i83mpqgrBag6ioCKG5OQGnMw6LZQ7Z2cvJyspGJ0znn/5wpMfrWjY1ltLSUqqqqnrUbg5Y63h8q5XbL1rMq3vKA+6iGijLr8voOzFeSSgFkMHUhvq76LW3t1NYWEhDQwMej4fS0lIqKiqoqqrC5XJhNptRSlFdXY3T6aSlpYXgYCOpqUbi4pwoVUJ8fDMRER0BFBw8tMEGXi9UVQVhtXYsHtfc3HF/T3LyEnJzl5OVlUV2dvYpMxS4XC7y8/P57YdHeG23lRsWp3Lzso7lsdPS0khISDilyS1QB4X4YxSWzJggznbjcqDDmagyj0W1fKB7JT45aeeFv33OVelGnnvzr3hK9uFxFPPB5wVEGNxkpycTGhpKU1MThwutOBpc5GUkEBcXidFYist1nJkzW4mPbyI8vJqwsBqMxqGtXNreDuXlHYMNKitDaWpKAKaQmbmKpUtXMHduItOnT+8RPpu25VPe7sFz6BBaa7TWHLDWcaKykZuXT8euothWG8XGm+fy0q4SrovJ6PcCHKiDQvzVtzDc+2WEONsEVE3pTHya7L3P7/zhIO8ctPnuSenaZjRC6qdv7+fp/9nB9XkJrMkwceDAAQ6cKGbrYRsX5KaTPSWRYnsj//vZcZZPjyHMHM6HhY1cPCsMS1AhoRF1eMKbmD+1lZjIaiIiGobc39PSoigrM1BeHtwZPvFERMxh7tx15OUtIDIykunTpxMSEuJ7jtfrpbi4mJqaGt/waq01B8vqefyjcjbdcxErZiYNew697o/1vmG567lj1Rw1HpdfFyIQjdvmuzNxw1r3Y7ywowjwz42SAB6Ph+LiYqxWK1arlRMnTlBeXk5FQwufFNQzPzOFLyqcrMuOJ8rQzolaNylxFmKCPTgcJ2nSFWRkKBLjm0hLqCMsvIZIy9D7e+rrO5rcbLYQ7HYTra3JJCYuJidnFTNnZvUZPi0tLb4bYaEjeLpmsE5PT+eNL+sGvXbNUC7m/hpcEujGw6SZQoyWcRtK0P+nydNd6IAhfaLtfoyuJpmhXCwqKirYs2cPNpsNq9VKVVWVbyaDmJgYYmJiMJlMVFZWYmts48+Hqlg5NQJPQxUuXY0zqJp1S6OYFFdHWJidqKg6QkNbh3yuqqo6mtxstmAcjnAgncTERWRnLyErK4uoqCgyMzN94aO1xuFwUFJS4mty6+rb6Vo0zmQy9Xms/mqyK2cm8Mf9ZaN2c+hEuHhLn5I4243LPiU4fd/CQPeMDPZ+kr6O0XsSRK01NpuNjz/+mLKyMmpra2loaKC5uaPWEhERwbRp04iKisJsNhMfH4/RaKS2traztuEkKspJTo4FU2w1P8h1kBTXQGRkPUbj0G4udbvBVmGgqMSAvcqEsyGK0NAZRETMZs6chSxblk10dDSZmZm+m0V/9f4xSprqmWpwcvDgQaBjhNvJSif3XrKA+fPnD3nRuL76RbpmExjNm0MnwkV7uPfLCHG2Caia0mA+TZ7uU/RgPmF33+fCVAtvbtvHQ795m7ZqG/NSLHxeZOectFASI0KJjIwkKyuLyZMnU1ZWRmlpKcHBwVRXV1NSUkJjYyNGo5ucnEimTTMSHFxOWJiD6Oh6zOYagoKGenNpx3xu5eVGiq0hHLeGUdmSxtyU2Vx43jl4Y9J48qMyNv39pZyfPQkAp9PpW+oB8M3n9mV5A0/sqgZAVcAAAAzXSURBVGfTN873+yfz3lMw+fvT/0SsKQlxthuXzXeD7Yc4XWdx98f+34UzqayspKSkhPz8fE6ePMlnBdXEW4Kx0Npxf48O4f3yIFYunMNlmSbe3/U57x4o57wpJlqqbTQ3N2MwGEhNjSQnJxKLpZqwMAeRkbWEhzsIDq4Z8ut0Og3YbMGUlRmpq4vE6YzDZMqivT2G3Nw8jEmZPLm9nLzsmXxZ2cSmDQvIjPBgtVp9tZ2rF6QAX9XY+lo0bjQu7t33+ZvthTy4bgZ3rfDf0gzSzCXExDQuQ2kwui5Ss5ItHCyr5z+uySY5yElVVRWb/7aH//30CAsmh3PAWs+aGVFMsoQQFRVFeno6U6dOpaqqivz8fOrq6qiqquKT/UcwBytiwkOxWCKYNSsBY0QjweG1LJzhxWgsIzjYhsHQMOSy1tYGU1ERgs0WSkNDNI2NsQQHT8frtZCVlcXy5cuJjY0lMzMTo9GI2+1my9a9PLJlD/+/vXuPrbK+4zj+/rZA77SFtlwL0gpaRJxSFTFMF4xTEnE6pxhRMGY6N5epG2HJ/tgys2TTzC3LMMiiYzNx081swzkvOG9TZgeEwcBwKVCx0ErvN4TS9rc/zqFW5NDnlHPO85xzPq+E5JQezvnyTft8z/P7fZ/vs/LLs5g7tYjthzp47JXd/HTZQm688sKox+kMPav5xW2fLmUmYoJ4JLHcG4zFe+oiU5H4S5mi1N8fOkM4dOgQb23by5MbdnLTBUV0dbbz4vYGwFhy0SQ6ejN48+Bx7lz0BRZOz+fd7Xt55vUtzC06wUBP++CGfn5+PhMnllJRkcOsWVlMmHCMrq4dZGQcYvTow2RkRHdxaX8/NDWPob4hi+amXHo6x9LdXUJm5gyaujNxueNZdc8tlJaWUlFRwahRo0ITrPvbuGBi/mCR2Vbfzv7mTyiZXM7l502NycEzVrfgOMmv22PEgs7AYkPFXUYq6YvSCy+8wK5duxgYGKCwsBAzY+P+Ns6dMp7qmVPp6+vjP7sOsvbFf5H7STNt7e2UF+cyuTiHvLw8xo4dS1bxOKZeMJHl107gwIE36e2tZfToBkaPbsQsuotLe3uNpqZsjjQXsHOf0dWez44Pc6goq2JScRG9Ywp4r2cSj9yxkFuunkdNXRv3Pfk6D11ZykVTiz7zWnvb+nj0vTZW33lpwq7HOjkj7SsXT+HtPU2+H4z92DfSXtXZU3GXkUr6orRhwwZaW1tpaWlh79691NbWDt6YLSMjg7y8PLKysmjNLKYufzrfu7WCS8fV0thYM2TJrQmz6P5/R49m0tKST3f3OA4fHs2RI9l0d48nM3MK2dk5FBcX01Y0k5cODPDgTVfynUUz2b9/P93d3Wyrb+exV3Zz/YUTeXlHI4/edRU3zL/gtEtusThARrsM9tBz/41p63Y08QzXlp+oC0l18erZU3GXkUjalvCT1q1bR2dnJ7m5uUybNo1Vq1Yxb94Mmpu3sGPH3+nr28cJt5uMzMMU5fUA0NMDBQXeXr+zc8zg7RPq6zM5ciSb9vZCsrMnkZWVTUFBAQsXLmT27NlUVlbS1dXFgQMHBidY3zp3Ak/97S3K+pu4+ap5FBQUUF0NHbnloYPeTdew5IrIB71YtD2fqUX+dF2Hb+9piusYn2hu8+DHWKGgjjJKNqnYsi/BEbii5Fw/x47V8eCD86mvfze83LaF3t4nqKkJXSOUk+PttQYGoKMjl66uYo4eLeWjj04Wn7Hk5JSQnZ0dLibVLFt2CRUVFTQ0NPDxxx9jZpgZnZ2dbN26lcLCQrpyJ7H6f808vfI2FlSW8LXwQXdm1XEWFBREddCLxQHS6zy1RN0bJmjx+P2eqUrFXeIpcMt3TU1/YefOm6P6NydOGB0dBbS0FdHQWsxATwENDWNoackhL6+YnJzQPtPcuXO57LLLmDFjBnV1dYMXwg5VXl5OWVnZaZfczrQ8Fe1Zy0jW5SO9/6/fqGXjvpaIS1KJ3pwebonMj81ybdDHhvaUZKSSdk+pp2cXmzZVnfZ7x45lcrAxl77jZdA/icbGrPAfo7Dw0+JTVVVFdXU15eXl1NfX8/ymg5xblsfcKaGGiaysLJqskD0tJ3yZPj7SA+TpDgj3PbMFgLsXnBOI9X3tN6Q2FXcZqaQtSgMDJ6ipqeTQIejqKqaxrZjXNh/DdeezdU83l8+awrSyYvLy8jj//POZM2cOkydPHryz6tBZbuPGjWP69OnU1LWlzKe7eA6TjVVsqZBnEYmtpC1KJ61atYrW1lZycnKoHyji3Ubjjqsv4o750z5TeDIyMgZvGncmsf4E7+cnxpPLY1dWjudbp2w0p9JtHkQkdSR9992KFSvo7u5me307r766m+XXTmdDvZ3xRnFnEuuOoWg6zWLp1E3mUy2oLPHtrOR0hcfPeEQk+UQ3JjqBqqqqODFuBqt3wG9XLuWxe2/gieXzeeDZrWzc1xz16516MB/Jaww1tNPs8dd2J2SZamjhe/ja8wbf/2z/LyIiQRHYogRnHvcfjXgdzIeefS27fFrczwhilY9orHl73+fytHFfM2ve3he39xSR9BXYPaVYitdeRzp0mql5QURiIekbHYIunQ7W6VB8RSS+vBalQC/fBZkfS2l+SfQypYikr8B23wVdOnWaaayMiCSKipJENT5JM+NEJJ60fCeDhedkl93JPaS5UwvTaplSRPynRgcB1MwgIvGlRgeJipoZRCQIVJQEiP3ECxGRkVBREo0vEpHASMmipNE40VEzg4gERUoWpTN1k8nnfeOqys/tIS2oLNHtJkQk4VLyOqWhE7zVTSYikjxS8kwJ1E0mIpKMUrYoqZtMRCT5pGRRUjeZiEhySsmipG4yEZHkpDFDIiISdxozJCIiScdTUTKz68xst5nVmtn34x2UiIikp2GLkpllAquB64HZwO1mNjvegYmISPrxcqZ0GVDrnNvvnOsF/gjcGN+wREQkHXkpSlOAj4Z8XR/+OxERkZiKWaODmd1rZpvNbHNTU1OsXlZERNKIl6J0CCgf8vXU8N99hnNurXOu2jlXXVpaGqv4REQkjXgpSpuAmWY2w8zGAEuB9fENS0RE0pGni2fNbDHwSyATeNo595Nhnt8EfHiWsZUAmgt0espNZMpNZMpNZMpNZLHKzXTn3LDLaHGZ6BALZrbZy9W/6Ui5iUy5iUy5iUy5iSzRudFEBxERCQwVJRERCYwgF6W1fgcQYMpNZMpNZMpNZMpNZAnNTWD3lEREJP0E+UxJRETSjO9FabgJ5GaWZWbPhb9fY2bnJD5Kf3jIzcNm9oGZbTezf5rZdD/i9IPXyfVm9lUzc2aWNp1VXnJjZreGf3Z2mtmziY7RLx5+p6aZ2ZtmtjX8e7XYjzj9YGZPm9kRM9sR4ftmZr8K5267mV0Sl0Ccc779IXTd0z6gAhgDbANmn/KcbwJrwo+XAs/5GXPAcvMlIDf8+H7l5nPPKwDeAd4Hqv2OOyi5AWYCW4Hi8NdlfscdoNysBe4PP54N1PkddwLz80XgEmBHhO8vBl4GDJgP1MQjDr/PlLxMIL8R+F348Z+BRWZmCYzRL8Pmxjn3pnPuaPjL9wmNgEoHXifXPwL8DDiWyOB85iU3XwdWO+faAJxzRxIco1+85MYBY8OPC4HDCYzPV865d4DWMzzlRuD3LuR9oMjMJsU6Dr+LkpcJ5IPPcc71AR3A+IRE569op7PfQ+hTTDoYNjfhpYVy59xLiQwsALz83MwCZpnZe2b2vpldl7Do/OUlNz8ClplZPfAP4NuJCS0pJOSOEaNi/YKSeGa2DKgGrvI7liAwswzgcWCFz6EE1ShCS3hXEzq7fsfMLnTOtfsaVTDcDqxzzv3czK4AnjGzOc65Ab8DSxd+nyl5mUA++BwzG0XolLolIdH5y9N0djO7BvgBsMQ5dzxBsfltuNwUAHOAt8ysjtD69/o0aXbw8nNTD6x3zp1wzh0A9hAqUqnOS27uAZ4HcM79G8gmNPtNPB6TzpbfRcnLBPL1wPLw41uAN1x41y3FDZsbM7sYeJJQQUqXfQEYJjfOuQ7nXIlz7hzn3DmE9tuWOOc2+xNuQnn5nforobMkzKyE0HLe/kQG6RMvuTkILAIwsypCRUk3iAtZD9wV7sKbD3Q45xpi/Sa+Lt855/rM7AHgVT6dQL7TzH4MbHbOrQeeInQKXUtoE26pfxEnjsfcPAbkA38K934cdM4t8S3oBPGYm7TkMTevAtea2QdAP7DSOZfyqw8ec/Nd4Ddm9hChpocVafIhGDP7A6EPKyXhPbUfAqMBnHNrCO2xLQZqgaPA3XGJI03yLSIiScDv5TsREZFBKkoiIhIYKkoiIhIYKkoiIhIYKkoiIhIYKkoiIhIYKkoiIhIYKkoiIhIY/wcpPX6GfXW1BQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with pm.Model() as model_robust:\n",
    "    family = pm.glm.families.StudentT()\n",
    "    pm.glm.GLM.from_formula('y ~ x', data, family=family)\n",
    "    trace_robust = pm.sample(2000, cores=2)\n",
    "\n",
    "plt.figure(figsize=(7, 5))\n",
    "plt.plot(x_out, y_out, 'x')\n",
    "pm.plot_posterior_predictive_glm(trace_robust,\n",
    "                                 label='posterior predictive regression lines')\n",
    "plt.plot(x, true_regression_line, \n",
    "         label='true regression line', lw=3., c='y')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There, much better! The outliers are barely influencing our estimation at all because our likelihood function assumes that outliers are much more probable than under the Normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- `PyMC3`'s `glm()` function allows you to pass in a `family` object that contains information about the likelihood.\n",
    " - By changing the likelihood from a Normal distribution to a Student T distribution -- which has more mass in the tails -- we can perform *Robust Regression*.\n",
    "\n",
    "The next post will be about logistic regression in PyMC3 and what the posterior and oatmeal have in common.\n",
    "\n",
    "*Extensions*: \n",
    "\n",
    " - The Student-T distribution has, besides the mean and variance, a third parameter called *degrees of freedom* that describes how much mass should be put into the tails. Here it is set to 1 which gives maximum mass to the tails (setting this to infinity results in a Normal distribution!). One could easily place a prior on this rather than fixing it which I leave as an exercise for the reader ;).\n",
    " - T distributions can be used as priors as well. I will show this in a future post on hierarchical GLMs.\n",
    " - How do we test if our data is normal or violates that assumption in an important way? Check out this [great blog post](http://allendowney.blogspot.com/2013/08/are-my-data-normal.html) by Allen Downey. \n",
    "\n",
    "Author: [Thomas Wiecki](https://twitter.com/twiecki)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": false,
   "threshold": "3",
   "toc_cell": true,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
